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Polynomial Approximations to Finitely 
Oscillating Functions 


By William J. Kammerer 


1. Introduction. Chandler Davis [1] established the following theorem: I[f 
Yo, ¥1, °** » Un are real numbers such that vo > v1, 0; < V2, V2 > Us, --- , then there 
bexists a unique polynomial P of degree n and a set of points yo, i, --- , Yn such that 


(1) P(y;) v; i=O0,1,---,n 
(2) P’(y:) = 0 ; 2.--> n=l 
O=-yw<n<---<y=1. 


The main result of this paper is an algorithm for the calculation of this polynomial, 
which is first motivated by an independent proof for the existence of P. 

A function f is said to be finitely oscillating if it has at most a finite number of 
Telative extrema. The following mode of approximating a continuous finitely oscil- 
lating function f in the uniform norm so that the oscillations are preserved, is dis- 
cussed: first obtain a polynomial P of minimal degree which has the same variation 
as f and then obtain an increasing polynomial Q such that P(@) agrees with f at 
all its relative extrema. Two theorems are given in the last section, to show that 
this method of approximation is always possible. 


2. Proof of Theorem. Let D and Dy, denote the set of all n + 1 tuples 
X = (%,%,°°* , 2n) such that O = m <4 < --- <a, = landO = »% <2 < 
+++ < 2, S 1 respectively. Let @ denote the class of all polynomials p of degree n, 
for which there exists an element X ¢ D, such that p(z;) = »;(i = 1, 2,---,m). 
A polynomial p ¢ @ can be written in the following form: 

(3) p(x) = [vo] + [xo , ri](2 — 20) 
+ on: + [to, 11, a. » tae _— Xo) eS (x — Zo 1) 


where the bracket function is defined by 


[zJj = »; 


(4) [rin1, «ae » Tigtl — [x;, coat » Vitk i] 


[vi Lis, oy Lisel nad 
Vitk — 


Lemma 2.1 If X ¢ Do and n 2 1 then [|X] > 0 for n even and |X| < 0 for n odd. 
Proof. The proof proceeds by induction. The lemma can be shown to hold for 

n = 1 by direct computation. Assuming it true forn = k — 1, one has |x», 1 , 

%-.|) < 0 (>0) and [z,, 22, --- , | > 0 (<0) for k even (odd). The lemma is 

therefore true for n = k by (4). 
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Lemma 2.2. a) If n is even (odd), the function |X| assumes its minimum (mazi- 
mum) in D. 

b) Let | [X]| attain its minimum value for X ¢ D, at Y = (yw, , 
* 5 Yn). Then the polynomial p(y.) = vue(k = 0, 1, --- , n) is the desired poly- 
nomial P. 

Proof. By lemma 2.1, [X] is of constant sign in D. The function | |X] | approaches 
infinity as X approaches the boundary of D, and is continuous in the compact set 
MteSu, M+ € S %,°*+ »La1 + € S 2 for arbitrarily small « > 0. 

The following notation is introduced, to simplify the proof of part b: 


a(p,i) = min {z:z e I(p,i)} F 
#=0,1,°---,n 
B(p, i) = max {x:z € I(p, i)} 
where 
; p(x) = v, if v; > 9; | ) 
I(p, i) = «| ; | and tia <2 < Lins} 
p(x) S vif 45 < vis | | 
#=1,2,---,n-—1 I(p,0) = 0 I(p,n) = 1 and pe@. 


Let Y be a vector in D at which |[X]| attains its minimum and let p be the poly- 
nomial in ® which satisfies p(y;) = v;. If p does not satisfy (2), then there exists 
an integer k, such that a(p, k) ¥ B(p, k). Let x% = }{a(p, k) + B(p, k)} and con- 
sider the following two polynomials 


h(x) + [Vlge(2, Y) 


p(x) 


Pilz) = h(x) + [Yo, +++ Yeo, Te, Yess, *** » Ynlge(z, Y) 


where k(x) is a polynomial of degree n — 1, and 
(5) g(x, Y)(a — ye) = I (a — yi). 


By construction one has | p;(2)| < | p(a,)|. Investigation of the possible cases 
contradicts the hypothesis that Y is a relative extrema of [X]. 

It should be observed that the above theorem and proof are valid for poly- 
nomials of the form > i» a,f(x)* where f is a strictly increasing differentiable 
function on [0, 1]. The proof is identical except for notation. 


3. An Iterative Procedure. 

Step 1. Choose an arbitrary element X, = (a , 2, -:- , 2n') in D. 

Step 2. With one of the standard interpolating formulas construct the polynomial 
pi € &, such that p:(a,') = %,k = 0,1,---,n. 

Step 3. Determine the vector X2 = (a0 , a1, --+ , 2») in D such that p,'(z,") = 0 
fors = 1,2,---,n— 1. 

We now have a new element X; ¢ D. To obtain p: , repeat this process beginning 
with step 2, using X2 in place of X, and making the obvious change in subscripts. 
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Continuing this procedure, we obtain a recursive process for obtaining a sequence 
{pi}. 

THEOREM 3.1. The sequence |p, converges uniformly to P. 

Proof. The proof will proceed by a series of lemmas. 

Lemma 3.2. If pe@, then 


Ito, *** » Ze-1, a(p, kh), tear, °° y Xn] = [t0, +>, t+, B(p, hk), tess, -°* , nl 


Proof. By (3), p can be written in the following two forms 
(6) p(x) 
(7) p(x) 


h(x) + [x0 , *** 5 Tk, a(p, k), Tet+i,y °°" 5 Tn] x(x, X) 


h(x) + [xo, “°° 5 Tht, B(p, k), Tktis*** ys Zn] 9u( x, X) 


where h(x) is a polynomial of degree n — 1 and g,(2z, X) is defined as in (5). To 
obtain the desired result, subtract (6) from (7). 

Lemma 3.3. Let p be an element of @ and let N(p) = maz;|m; — v;| where 
m; = maz | p(x) | for xe la(p, i), B(p, t)| and i = 1, 2,---,n — 1. Then the 
following inequalities hold for the sequence {p,j : 

a) If N(pi) # O then | [Xiu] | < | (Xd | 

b) | pix) — piss (x) | S | (Xa) — [Xia] | for every xe \0, 1) 

c) N(pi) S mazogesi| p(x) — piss (x) |. 

Proof. Let i be any positive integer and let {hy}, k = 1, 2,---,n — 1 be the 
polynomials of degree n such that 





hy(tn) = Um form = 0,1, --- ,k 
hi(Bm) = Um form =k+1,---,n 


where 8», = 8(p;, m). Let Z, denote the vector (2 Ti? Wey sony »*** , Bn)- 
Then J(hy ,m) C I(hy4,, m) form = k + 2, --- , n, for otherwise hy(x2) — hyss(x) 
would be identically equal to zero, instead of having n — 1 simple roots in {0, 1). 
As in the proof of lemma 2.2, we have | [Zo| | = | [Z:)| 2 --- = | [Z,] | with at 
least one of the inequalities being strict, since by assumption, N(p,;) # 0. This 
proves part a. 

For k = 0, 1,---, m — 1 one has| hy (2%) — Pyar (2) | = | [Ze] — [Zesal| - 
ge (a, Ze) | S | [Ze] — [Zess] | . Now apply the triangle law to obtain b. 

The proof of part c follows from the existence of an integer k, such that N(p,;) = 
| pila’) — Diva (a,'*") . 

To complete the proof of theorem 3.1, observe that lemmas 2.1 and 3.3.a imply 
that {[X,]} is a Cauchy sequence. Lemma 3.3.b implies {p,} also forms a Cauchy 
sequence in the uniform norm on [0, 1] and therefore convergent to a polynomial 
P of degree n. Part c of lemma 3.3 implies that P satisfies conditions 1 and 2. 

This iteration can be carried out by the use of standard subroutines available in 
most computer libraries, and from all empirical evidence the convergence seems 
quite rapid. To illustrate, we shall calculate the third-degree Chebyshev poly- 
nomial on [0, 1] by the use of this method, starting with X, = (0, .25, .5, 1). 
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pi(x) = 1 — 222 + 682° — 482° 

X2 = (0, .20723911, .73720533, 1) 
p(x) = 1 — 18.4005412 + 48.9709052" — 32.5703642° 

X; = (0, .24979357, .75256767, 1) 
p(x) = 1 — 18.000432x + 48.0021382” — 32.0017072° 


X4 = (0, .24999339, .74999783, 1) 
p(x) = 1 — 17.9999992 + 47.9999992" — 31.9999992" 
4. Approximation by Composition. 

THEOREM 4.1. Let X = (2, %1,°**,2n), ¥ = (Yo, fi, °** 5 Yn) be any two ele- 
ments in D. Then there exists a polynomial Q such that Q’(x) = 0 on [0, 1] and 
Q(2z;) sa Yi, t - 0, 1, ee ee 

Proof. Define the elements Z; and Z, in D as follows: 


Zi = (20, 215 °** 520) = (Yo, Hye — m1) +, HYys — yr) + ye, -** 5 Yn) 
Z2 =“ (zo, a’, << = Zn) = (yo, i(% hej Yo) = is i (Ye sal 1) = 2s, *** » Ba) 


Let f; for i = 1,2,3,--- , 2” be distinct piecewise linear functions, which are linear 
on the intervals [y;. , y;] i = 1, 2, --- , n and such that f(a.) = z' or z%, k = 0, 
1, --:, n. Define « = min; min, |2z/ — y;| fori = 1, 2,---,n — 1,k = 1,2. 


Using the fact that the Bernstein polynomials of a continuous increasing function 
are increasing and uniformly convergent, there exist increasing polynomials Q;(x) 
on [0, 1] such that | f(x) — Q(x) | < }e fori = 1, 2,---,2”"” (see Lorentz [7] 
p. 20-23). The vector Y is contained in the convex hull of the vectors (Q;(2»), 
Qi(ai), --- , Qi(an)), i = 1,2, --- , 2”, and therefore there exists a convex linear 
combination of the Q,’s which will give rise to a desired polynomial Q. 

THEOREM 4.2. Let f be a continuous finitely oscillating function on \0, 1| and let 
€ > 0 be given. Then there exist polynomials P(y) and Q(x), such that 

a) f(x) and P(y) are equal at their corresponding relative extrema. At the relative 
extrema of f, 


. " adPli)) _ 
P(Q) =f and — 0. 


b) The polynomial Q is increasing and | f(x) — P(Q(x)) | < eon (0, 1). 

Proof. Let the partition X = (a, %1,+-+, 2») be the points of the relative 
extrema of f on [0, 1]. By a previous theorem, there exists a polynomial P(y) and a 
partition Y = (yo, ¥1,°*:, Yn) such that P(y;) = f(x), i = 0,1,---, m and 
P’(yi) = 0,7 = 1,2, ---,n — 1. Let X’ be a refinement of the partition X which 
satisfies the following condition 2; = ta < t2 <-++: < 2iw; = %i4. such that 
Ni|f(ris) — f(xisns) | = |f(as) — Slain) | S Niefor i = 0, 1,---,n — 1, 
j = 1,2,---,N;. Define Y’ to be the refinement of the partition Y such that 
Yis © (Yi, Yess) and f(2;;) = P(yi;) fori = 0,1,---,n—1,7 = 1,2,---,N;. By 
theorem 4.1 there exists an increasing polynomial Q() on [0, 1], such that Q(2,;) = 
y:;, and therefore by construction | p(Q(x)) — f(x) | < eon {0, 1). 

As a concluding remark, let C represent the class of all composite polynomials of 
the form P(Q), where both P and Q are of degree greater than unity and Q’ = 0 
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on [0, 1]. Not every polynomial of degree four or larger can be so written (see 
H. Levi [6]). However, since every continuous function on [0, 1] can be uniformly 
approximated by a polynomial, (i.e., a finitely oscillating function) one finds that 
the completion of C in the uniform norm on [0, 1] is the set of all continuous fune- 
tions on [0, 1]. 

The question of the existence of variation preserving approximations arose from 
the investigation of syntonic functions by Professor P. C. Hammer in [2] and [3}. 
I would also like to thank the referee for pointing out that the ideas in references 
[4], [5] and [8] bear a relation to the problem treated here. 


Georgia Institute of Technology 
Atlanta, Georgia 
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Table of a Weierstrass Continuous 
Non-Differentiable Function 


By Herbert E. Salzer and Norman Levine 


Many studies have been made of continuous non-differentiable functions [1], the 
most famous of which is Weierstrass’s W (a, b, x) defined by 


(1) W(a, b, z) = >> a” cos (b"x2), 0 <a < 1, b an odd integer. 
n=l 

It is shown in some books [1], [2] that for 

(2) ab>1+ a 


W (a, b, x) is continuous everywhere and has no derivative anywhere, but Bromwich 
[3] improved this condition to 


(3) ab >1+ 5% (1 -a), 

which, according to Hardy [4] is the sharpest result (as of 1916) for no derivative, 
finite or infinite. (Hardy showed b > 1, ab 2 1 sufficient to establish the non-exist- 
ence of any finite derivative. He also showed that those same conditions, together 
with a(b + 1) < 2 for b = 4k + 1, permitted the existence of an infinite derivative 
at certain points.) To illustrate the difference between (2) and (3) for a = 3}, 
(2) requires b 2 13, while (3) permits b = 7. However, as far as the authors know 
there may be considerable work to be done in the direction of lowering the bound 


of 1 + os (1 — a) in (8) for the case of no derivative, finite or infinite. 


Owing to the unusual nature of W(a, b, x) and the absence of any previous table, 
or even graph, despite the countless number of theoretical papers, it was believed 
that an extensive table of this Weierstrass function for some typical pair of param- 
eters a and b might be of value as more than a mere curiosity, namely for suggest- 
ing or motivating further research, and for its interest to workers in numerical 
analysis. Thus, in this last connection, it might be of interest to determine empiri- 
cally what results in numerical integration and possibly interpolation are available 
from the continuity alone. That W (a, b, x) is integrable follows from its continuity, 
and one might be curious to see the results of applying standard numerical integra- 
tion formulas where the usual derivative formulas for the remainder would be 
inapplicable. Likewise, one might be curious to test out standard Lagrangian inter- 
polation, where the remainder is often expressed in terms of derivatives. (We can 
write down interpolation and numerical integration formulas, avoiding derivatives 
in the remainder terms by employing divided differences and integrals with divided 
differences in the integrand, respectively. However, one usually estimates divided 
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glancing at the results of standard numerical differentiation and interprétation of 
the results in the light of the knowledge that W(a, b, x) has no derivative. 
For tabulation of any W(a, b, x), it is immediately apparent from (1) that 


(4) W(a,b,1 +2) = — W(a,b, 2), 
so that the range of x need not go outside (0, 1). From (1), 


W(a, b,0) = —W(a,b,1) = a/(1 — a); 


(5) 
W(a, b, 3) = 0. 
From the trigonometric identity 
(6) cos(mr(3 + t)) = #(—1)°”” sin mat, m odd, 
we have 
(7) W(a, b,} +t) = —W(a, b, } — 2), 


so that for complete tabulation of any W(a, b, x) it suffices for x to range from 
0 to 4. 

In connection with the choice of a and b, it is apparent that for a close to 1, we 
can choose b as low as 3, but the convergence of the series in (1) would be too slow 
for practical calculation of W(a, b, x) to high accuracy. Making a very small would 
give rapid convergence, but for accuracy fixed at a certain number of decimal 
places as a tends to get very small, say 


a=.6, w>n=h+a-o} fe 


becomes enormous and W(e«, b, x) becomes essentially the first term of (1), € cos 
(b"xx), whose graph would appear like that of a very highly oscillatory function of 
small amplitude. As a compromise between these two extreme types, we took 
a = 3 and b = 7. The choice a = 3 did not lead to too many terms of (1), 50 terms 
giving a truncating error < 4-10°"°, and yet there were sufficient terms beyond the 
first few to give a graph that is characteristic of W(a, b, x) rather than a predomi- 
nantly sinusoidal type of curve. The b = 7 barely satisfies (3), thus tending to 
minimize the oscillatory behavior of W (a, b, x) and to facilitate graphing. We shall 
denote W (a, b, x) which is tabulated here for a = 3 and b = 7 by W(z). 

This present table of W(x), « = 0(.001)1 to 12D, was printed out and rounded 
from a preliminary calculation on the IBM 704 to several more places. Two sepa- 
rate and independent print-outs, supposedly identical, were proofread against each 
other, with just a single print-out error turning up. Naturally, no differencing check 
could be made upon the correctness of this table of W(x), but every value under- 
went the following final functional check: 


(8) W(7x) = 2W(x) — cos (7272), 


which was performed by desk calculation upon W(x) on one of the preliminary 
print-outs. The results showed W (x) to be correct to around 14D. In employing (8), 
W (72x) was found in the table as + W(2’) for some suitable zx’,0 S x’ <S 3, accord- 


ing to (4) and (7), and cos (7zx), after reduction of 7x to the first quadrant, was 
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Fig. 1.—Illustration of a Weierstrass, Everywhere-Continuous Nowhere-Differentiable 
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looked up in a well-known 15-place table at intervals of 0.01° [5]. The final 12-deci- 
mal table was checked by reading it several times against one of the print-outs, 
and it is believed to be correct to well within a unit in the 12th decimal. 

The purpose of the accompanying figure, which is merely a broken line graph of 
the table of W(x), is to furnish at a glance a view of the peculiar behavior of W(z). 
Of course, the graphical picture would be more complete if the time and means were 
available for calculating W(a, b, x) as a function of a also, and for a sequence of 
permissible odd integral values of b (according to (3)) to correspond to each a. 
Although no offhand justification could be found for drawing anything smoother 
than a broken line connecting these 500 points, one still finds its ripples of irregu- 
larity, superposed upon a broader pattern of smoothness, to be quite revealing as 
to the nature of W(x) and how it might appear under repeated “magnification” 
(i.e., subtabulation). 

To establish (8), replace x by 7z, in W(x) = }-%_, cos (7"xx)/2", to get 
W(7x) = 2 > cos (7"*'xx)/2"*" = 2 ay cos (7" xx)/2” = 2W(x) — cos(7xz). 
By repeated application of (8), 

W(7"x) = 2W(7" x) — cos (7"xx) = 4W(7" x) — 2 cos (7" xx) — cos (7"x2r) 


= 8W(7" *x) — 4cos(7" xx) — --- ete. until we reach 
(9) W(7"x) = 2"°W(z) - > 2 cos (7"” xz). 
From (9), for = 1/7", W(1) = — 1 = 2°W(1/7") — Fe cos (2/7"), from 
r=0 
which 
n—l 
(10) W(1/7") = — 1/2" + > cos (x/7’)/2””. 


r=1 


Letting n — © in (10), we see at once that 


(11) lim ‘= 2’ cos (x/7)} [2 = 1 


n>o r=1 


To test the value of standard numerical integration formulas upon W(x), whose 
integral is given by 


(12) [ wit) dt =+ 1 > sin (7" xx) /14", 


0.1 0.2 . 
the values of [ W°(t)dt, | W (t)dt,---, f= W (t)dt were computed analyticaily 
0 0.1 0.4 


from (12), and then were computed numerically by both trapezoidal and Simpson’s 
rules at intervals of 0.001, with the following results: 





Interval True Value | Trapezoidal Rule ~ Deviation Simpson’s Rule Deviation 
| 





| 0.01899 29 0.01898 76 | —0.00000 53 0.01901 44 | +0.00002 15 
| —0.04145 65 | —0.04143 80 | +0.00001 85 | —0.04145 43 | +0.00000 22 
0.03084 62 | 0.03084 43 | —0.00000 19 | 0.03085 14 +0.00000 52 
0.00337 70 | 0.00342 54 | +0.00004 84 | 0.00340 27 | +0.00002 57 
—0.03298 02 | —0.03300 67 


—0.00002 65 | —0.03288 27 | +0.00009 75 
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The results show no recognizable advantage in Simpson’s rule. In fact, the sum of 
the absolute values of the above deviations in the trapezoidal rule is around 10“, 
while the sum of the absolute values of the Simpson deviations is around 1}-10™~. 
This may indicate that no higher-point formula will improve over the trapezoidal 
formula. 

Lagrangian polynomial interpolation at intervals of 0.002 was tried for the 2- 
through 7-point cases, for a mid-interval (i.e., already tabulated) value of W(x) at 
two different places, = 0.007 and x = 0.037, where the true value to 5D is 0.60807 
and 0.43362 respectively. At each place the error in almost all cases ranged from 
around 0.01 to 0.05. More specifically, for « = 0.007 the error fluctuated between 
0.01 for every even-point interpolation and 0.014 to 0.049 for various odd-point 
interpolations, and for z = 0.037 there were deviations of 0.032 and 0.055 for respec- 
tive 2-point and 3-point interpolation and deviations ranging from 0.001 to 0.021 
in the higher-point interpolation. On the basis of these two tests alone it would 
appear that one could not really count upon any systematic improvement beyond 
linear interpolation. 

Finally, out of pure curiosity, 2- through 7-point Lagrangian differentiation, for 
the “first derivative,” was tried out at the tabular interval of 0.001, for x = 0.002, 
and surprisingly enough, outside of the 2-point answer of —74 and the 3-point 
answer of —133, the remaining four cases all came within 6 units of — 150. 

From a casual look at the graph of W(z), it is apparent that in place of the de- 
rivative there is a general directional trend from any point 2 if we do not go too 
far away from 2, and we might seek a suitable quantitative estimate for an “aver- 
age slope” between 2 and x + h. (The discussion here is concerned with a suitable 
generalization of the left- or right-hand derivative, rather than the derivative.) 
One suggestion that would appear natural for W(z, a, b), or any other continuous 
function, would be to investigate the possibilities of the average of the difference 
quotient {f(z) — f(xo)}/(% — 2), which exists and is itself continuous for every 
x except 2% in the open interval (2, x + h). This average difference quotient or 
Dif(x%o) might have the following definition (assuming that it exists in the first 
place) : 


1 zoth 
(13) Dafa) = Ff” {Lf(e) — Hlae)I/(@ — 2)} de. 
That (13) may be a suitable generalization follows from the fact that when f’(2») 
exists, (13) exists, and 
(14) lim Dif (20) = f’(x0). 
This is seen at once from the replacement of {f(2) — f(xo)}/(x — a) by f’(ae) + 


¢ (x) in (13) and the continuity of ¢ (x) in the closed set (xo , x» + h) which makes 
e (x) integrable. Thus (13) exists and 


11 zoth 
Fal e(z)dz|- 0 as h-0, 
2) 





which implies (14). 

It is not difficult to find examples of continuous functions f(z) where f’(xo) does 
not exist and (a) also D,f(zo) does not exist, or (b) Daf(xo) exists but lims.o 
Daf(x%o) does not exist. But we may also have (c) no f’(2), with both Dsf(xo) 
and lims.o Daf(xzo) existing. In other words the existence of limyj.o Daf(zo) still 
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does not imply the existence of f’(zo). Such a counter-example,* whielr is due to 
the referee, is the following. Let z = 0 and 
f(z) =2 sin = (x # 0) 
f(0) = 0. 
This continuous function has no derivative at z = 0, but 
lim D,f(0) = 
h>0 


m=} [sin (2) ae 


exists since the integrand is bounded and continuous except at one point. This 
suffices. To estimate D, we let 


l/ne : 1 (n+1)4 1 
a i = (2) ae vr “4 puny dy. 


By the mean value theorem 


First 


= (—1)"-2/0,’ 
where 
ner <0, < (n+ 1)z. 
Suppose that h = 1/(n + a)x,0 S a < 1. Then 


(n+l) 
= (n+ ade | y* sin y dy + Toma + Inve + -~ |, 


n+a)r 
and therefore | D, | < (n + a)x|I,| <2(n + a)x/n'x° 


Therefore as h — 0, Dp also — 0. 
The authors wish to acknowledge the assistance of Mrs. Charlene M. Janos in 
checking the entire table of W(x) by the functional check (8). 


Convair-Astronautics 
San Diego, California 
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* Another counter-example found after that of the referee is the following: f(z) = 
xo(x), x ~ 0, f(0) = 0, where ¢(z) = 1 except in the intervals [(1/n — 1/n*), 1/n), within which 
¢(z) = 0. Now f(z) is continuous at z = 0 and has no derivative there. But 1/h fe ¢(z)dz — 1 as 
h— 0, because the “dipped-out” area becomes an infinitesimal fraction of the whole (also in- 
finitesimal) area between 0 and h, since as h ~ 1/n, we remove » poe 1/m? ~ 1/2n? ~ O(A). 
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TABLE oF W(x) = 3 cos (7"x)/2”" 








W 


(x) 


x 


W (x) 





.041 
.042 
-043 
.044 


.046 
.047 
.048 
.049 





1.00000 00000 00 


.80391 58298 49 
.61188 60438 58 


.53777 
-64747 
.87163 
. 76687 
- 60807 
- 43502 
-40541 
-56641 
-54275 
- 50694 
.34801 
. 22473 
.27196 
- 25665 
. 34500 
. 29744 
. 19896 
. 16232 
.07772 
- 20584 
- 28363 
31741 
. 28730 
. 11054 
.17279 
. 29881 
-48372 
54441 
. 32388 
. 26990 
33225 
.58370 
-74931 
56632 
.43361 
. 36496 
55435 
- 75565 
.66141 
.54244 
36553 
41175 
.54502 
.52178 
.49248 
. 30088 
22797 


60375 
48039 
69853 
71957 
34552 
94075 
06494 
31472 
36720 
91215 
25245 
91530 
45026 
87904 
48434 
09740 
02842 
54753 
75335 
44795 
56796 
47365 
16038 
29341 
94307 
59987 
84610 
21945 
99122 
13283 
74462 
92580 
30151 
09611 
86486 
26383 
98106 
42269 
61182 
67604 
53065 
90597 
46829 
60105 
88676 
64437 
19914 


27 
38 
23 
75 
61 
78 
76 
93 
27 
98 
87 
39 
68 
18 
78 
20 
26 
01 
97 
34 
36 
60 
97 
42 
92 
33 
58 
09 
78 
84 
04 
29 
91 
85 
58 
50 
40 
37 
70 
80 
06 
74 
96 
03 
02 
50 
12 


1.000 
.999 
-998 


.996 
-995 
-994 
-993 
.992 
.991 
-990 
- 989 
-988 
- 987 
-985 
-985 
- 984 
983 
. 982 
-981 
-980 
.979 
.978 
977 
.976 
-975 
.974 
.973 
.972 
971 
-970 
-969 
-968 
. 967 
-966 
.965 
. 964 
.963 
-962 
-961 


-959 
.958 
957 
956 
-955 
.954 
953 
-952 
951 


-050 
.051 
-052 
053 
-054 
055 
056 
057 
-058 
059 


-061 
-062 
-063 
064 
-065 
-066 
067 
.068 
069 
.070 
O71 
.072 
.073 
.074 
-075 
.076 
.077 
.078 
.079 
-080 
-081 
-082 
-083 
-084 
-085 
086 
.087 
.088 
089 


.091 
.092 
-093 
.094 
-095 
.096 
-097 
-098 
-099 








. 23088 
. 20682 
. 27128 
. 16118 
-08069 
-02066 
. 11450 
.02295 
-01951 
.01151 
.09698 
.27187 
. 23653 
. 16965 
.01498 
.00239 
.20181 
. 25856 
. 21932 
.05550 
. 15690 
.01436 
.09812 
. 15074 
-09240 
. 24890 
- 20632 
. 11462 
.02304 
-09557 
. 19794 
. 22531 
. 20876 
-05397 
-04851 
.02128 
-03684 
.08433 
.01214 
.05215 
. 19576 
. 26338 
-21893 
. 22931 
. 19289 
. 36048 
.51624 
.54350 
-50350 
. 33848 


91433 ; 


52628 
71570 
34941 
56769 
12990 
55193 
65257 
72464 
68818 
14952 
35472 
00063 
63244 
43634 
70885 
95236 
31395 
57817 
33815 
95326 
83992 
03304 
56668 
88499 
58340 
59257 
33882 
18919 
81653 
01773 
98834 
59176 
57757 
66043 
19742 
58507 
37682 
33731 
21171 
22844 
74226 
60178 
50089 
54543 
94459 
68304 
09214 
10354 
06786 


39 
31 
71 
70 
04 
73 
19 
55 


72 
45 
21 
87 


74 
23 
04 


47 
20 
88 


35 
07 
00 
35 

















49 
19 
20 
94 
43 
63 
78 
18 
0 | 
42 

01 
41 
81 





58 | 
34 

76 | 
42 | 
87 

40 | 
76 | : 


.950 
.949 
.948 
.947 
.946 
.945 
.944 
.943 
.942 
.941 
.940 
.939 
.938 
.937 
.936 
.935 
.934 
.933 
.932 
.931 
.930 
.929 
.928 
.927 
.926 
.925 
.924 
.923 
.922 
.921 
.920 
.919 
.918 
.917 
.916 
.915 
.914 
.913 
.912 
911 
.910 
.909 
.908 
.907 
.906 
.905 
.904 
.903 
.902 


901 











—W (x) 








—W (x) 
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TABLE or W(x)—Continued 





















































a ee eS Oe 














x W (x) | = W(x) 
100 | —.42532 54041 76 | .900 .150 | —.60928 87419 74 .850 
101 | —.60122 07728 99 | .899 .151 | —.48052 65499 66 | .849 
102 | —.71436 04664 23 | .898 .152 | —.49741 97079 36 | .848 
103 | —.69032 26595 69 | .897 .153 | —.49479 75145 65 | .847 
104 | —.43794 73064 53 | .896 .154 | —.47354 23522 65 | .846 
105 | —.40215 63534 91 | .895 .155 | —.49291 96963 68 | .845 
"106 | —.50671 30937 48| .894 | .156 | —.36979 28855 91| .844 
.107 | —.65237 17461 67 | .893 .157 | —.30677 09741 35 | .843 
108 | —.68741 46475 51 | 892 .158 | —.21917 49907 16 | .842 
.109 | —.44815 12393 09 | .891 .159 | —.21531 69983 41 | .841 
110 | —.33948 90492 28 .890 .160 | — .32774 87639 11 | 840 
111 | —.32696 83696 91 | .889 .161 | —.30751 75250 62 | .839 
.112 | —.42541 62768 18 | .888 .162 | —.27020 56659 27| .838 
113. | —.50701 25026 15 | .887 .163 | —.11494 78174 84 .837 
114 | —.36527 04053 46 | 886 .164 | —.05693 76441 49 | .836 
.115 | —.28717 02983 14 | .885 .165 | —.19542 29668 60 | .835 
.116 | —.19435 86539 72 884 .166 | —.30152 83944 41 | .834 
117. | —.20343 01549 89 | .883 .167 | —.38366 39179 31 | .833 
118 | —.27648 95287 $6 | .882 .168 | —.22909 39440 73 | .832 
119 | —.24091 87061 61 | .881 .169 | —.09500 91416 95 | .831 
120 | —.27427 89580 66 | .880 .170 | —.17303- 80098 99 | .830 
.121 | — .19594 53705 68 | .879 .171 | —.33806 67973 86 | .829 
.122 | —.14745 62719 19 | .878 .172 | —.55064 02832 16 | .828 
.123 | —.16077 90760 02 | .877 .173 | —.46584 78110 99 | .827 
.124 | —.16313 11716 57 | .876 .174 | —.30129 83741 07 | .826 
125 | —.30795 98441 70| .875 .175 | —.27803 69565 60 | .825 
126 | —.32779 07918 31 | .874 .176 | —.39280 65936 50 | .824 
127 | —.30642 17906 98 | .873 .177 | —.65182 79852 26 | .823 
128 | —.25457 54992 71 | .872 .178 | —.65182 39388 20| .822 
"129 | —.20521 28645 40 | .871 .179 | —.53351 02371 89 | .821 
.130 | — .38226 57006 18 | .870 .180 | —.44268 44274 72| .820 
131 | —.51008 60535 12 | .869 .181 | —.43578 01021 83 .819 
.132 | —.58897 88293 07 | .868 .182 | —.61985 87140 15 | .818 
.133 | — .51605 03218 22 | .867 .183 | —.64922 21692 63 .817 
.134 | —.37228 02555 49 | .866 .184 | —.62277 87433 01 | .816 
135 | —.48222 40135 76 | .865 .185 | —.54441 65431 73 | .815 
.136 | — .64476 90942 44 .864 .186 | —.43416 51101 85 | .814 
.137 | —.82656 60891 90 | .863 .187 | —.47153 43722 97 | .813 
.138 | — .78900 50242 13 | .862 .188 | —.44187 85912 52| .812 
.139 | —.58460 20582 18 | .861 .189 | —.48324 34653 21 | .811 
140 | —.58017 61018 65 | .860 190 | —.48100 45544 76 | .810 
141 | —.67358 93294 65 | .859 .191 | — .36210 24639 76 | .809 
142 | — .88334 97740 78 | .858 .192 | —.28385 07250 58 | .808 
143 | — .90195 54475 25 | .857 .193 | —.13957 34378 96 | .807 
144 | —.71735 67984 31 | 856 .194 | —.17132 73511 09 | .806 
.145 | —.63542 71888 15 | .855 .195 | —.24380 91207 97 | .805 
146 | —.60172 84929 47 | .854 .196 | —.21870 72532 52 | .804 
.147 | —.73979 05811 32 | .853 |- .197 | —.13838 71322 99 .803 
.148 | —.77996 61358 53 | .852 .198 |  .09801 82360 71 | .802 
149 | —.67821 23623 64 | .851 .199 | .14666 52327 61 | .801 
| —W (x) x —W (x) x 
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TABLE or W(x)—Continued 


























x W (x) x W (x) 

. 200 .06366 10018 75 .800 . 250 .70710 67811 87 .750 
.201 — .04175 22364 04 .799 .251 .59986 16383 91 .749 
. 202 — .07476 70750 18 .798 . 252 .42999 19525 7i .748 
. 203 . 16401 55220 04 .797 .253 . 36660 93235 94 .747 
. 204 .29971 35815 61 . 796 . 254 .32795 62544 06 . 746 
.205 . 28236 84375 00 .795 . 255 .45218 11134 91 .745 
. 206 .10035 12674 93 .794 .256 .43000 78607 24 .744 
.207 — .08006 36504 29 .793 .257 .38270 70934 21 .743 
. 208 .06602 53228 00 . 792 .258 .32431 23791 08 . 742 
.209 . 22194 83909 79 .791 .259 . 18040 34507 72 .741 
.210 .29866 93766 43 .790 . 260 .20082 17490 15 .740 
.211 . 14347 00102 07 . 789 .261 . 19502 92282 86 .739 
.212 — .11156 27605 97 . 788 - 262 .28277 11533 62 .738 
.213 — .08927 61047 33 . 787 . 263 .33124 94143 10 .737 
.214 — .00503 43683 93 . 786 . 264 .19535 83704 39 . 736 
.215 . 12393 18196 44 . 785 . 265 . 13130 01934 73 .735 
.216 .07237 02520 68 . 784 . 266 .05923 97001 64 . 734 
.217 — .12870 66881 32 . 783 . 267 -20320 81841 97 .733 
.218 — .17190 63119 85 . 782 . 268 .38097 30175 89 .732 
.219 — .20539 46224 39 .781 . 269 .36488 72388 22 .731 
. 220 — .10396 93551 29 .780 .270 .30069 58598 63 .730 
.221 — .05670 84900 53 .779 .271 . 12888 96220 58 .729 
. 222 — .11058 04979 98 .778 .272 .21930 88571 52 .728 
.223 — .11053 47920 41 are .273 -45450 66523 44 .127 
. 224 — .22263 98542 44 .776 .274 .58788 37027 28 .726 
.225 — .20433 20524 80 .775 .275 .61063 78772 02 .725 
. 226 — .13899 14773 88 .774 .276 .37897 11709 32 | .724 
.227 — .05338 52274 30 .773 .277 .35491 00906 37 . 723 
. 228 .07179 91975 44 44a .278 .53317 01360 59 . 722 
.229 — .01896 62898 75 .771 .279 .74080 38401 87 .721 
. 230 — .07113 29711 74 -770 . 280 .87388 44641 26 .720 
.231 — .08351 49064 96 . 769 .281 .66344 41290 97 .719 
. 232 .03967 20418 75 . 768 . 282 .55360 44310 88 | .718 
. 233 . 28470 23445 34 . 767 . 283 .59858 96747 52 | .717 
. 234 .30373 06421 48 . 766 . 284 .75311 92971 19 | .716 
- 230 .25331 56946 85 .765 . 285 .93575 68089 02 | .715 
. 236 . 12247 07877 84 . 764 . 286 .80593 31523 57 | 714 
.237 .16091 46178 36 . 763 . 287 .70250 54785 63 | .713 
.238 .43416 31698 35 . 762 . 288 .62769 78735 12 | .712 
. 239 .57254 03757 22 .761 .289 .64051 31524 34 .711 
. 240 .60020 74057 93 .760 .290 .76998 70795 56 | .710 
.241 .39293 37880 48 .759 .291 .71343 52538 50 | .709 
. 242 .29091 21091 20 .758 . 292 -70111 02427 36 | .708 
. 243 .47723 86339 28 .757 . 293 .59700 86384 17 | .707 
. 244 .65452 02702 15 . 756 . 294 .48196 39792 70 | .706 
.245 .78049 97326 19 155 .295 .48841 28610 21 | .705 
. 246 .58771 13946 18 . 754 . 296 .43820 32711 49 | .704 
.247 .39392 08421 96 .753 .297 .53246 89138 76 | .703 
.248 .43534 55892 80 . 752 .298 | .50028 47645 04 | 702 
. 249 .53704 70311 88 451 .299 | .36415 33185 00 
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TABLE or W(x)—Continued 






































2 th We) x Ws) 
.300 .26286 55560 60 | .700 | .350 .00778 77869 67 | .650 
.301 14582 98589 87, .699 | .351 | —.17204 22086 50| .649 
.302 .28563 97407 01 | .698 .352 | —.23029 21241 86| .648 
.303 .36633 50157 87 | .697 .353 | —.23647 30864 44| .647 
304 .33282 81740 47 || .696 | .354 | —.01762 25274 42 .646 
305 21458 96315 30 | .695 | .355 .06669 22195 13) .645 
306 00939 20926 45 | .694 .356 .00075 62690 21 | .644 
.307 .10710 80078 80 | .693 .357 | —.02063 75236 12| .643 
.308 .25624 91704 22| .692 .358 | —.08880 13434 64 | .642 
.309 .37838 97945 27| .691 .359 .03980 77688 40 | .641 
.310 34385 20085 52) .690 | .360 10006 60841 69 | .640 
311 .09875 62075 19| .689 | .361 .12219 89553 34) .639 
312 .10737 46483 70| .688 | .362 .18737 55558 83 | .638 
.313 .23160 40973 59 | .687 .363 10954 21317 00 | .637 
314 .45535 62002 52| .686 | .364 .12819 87330 70 | .636 
B15 .54102 65479 93 | .685 365 .07610 38829 89 | 635 
.316 .33736 28357 70 | .684 .366 09741 56111 37 | .634 
.317 .28355 86706 20) .683 367 |  .22835 66376 41 | .633 
.318 | .30968 52189 51 | .682 .368 19796 35878 29 | .632 
.319 |  .51427 37371 70| .681 .369 .16937 55518 30| .631 
.320 .66458 80166 07 | .680 | .370 | —.00562-33217 53| .630 
321 .55383 02222 36| .679 | .371 | —.07189 57360 01 | .629 
322 .51306 08367 00 | .678 | .372 .04604 08396 66 | .628 
.323 | .43864 14186 98 | .677 .373 .07904 58408 70 | .627 
324 | .52349 24425 12| .676 | .374 .09708 84764 28 | .626 
325 63004 29627 66 | .675 | .375 | —.12756 11441 22| .625 
.326 .59308 72965 01 | .674 .376 | —.30043 50117 72| .624 
.327 .62794 97613 07 | .673 377 | —.27256 67762 10| .623 
.328 51805 41271 03 | .672 878 | —.21256 78517 98 | .622 
329 | .47323 47350 89 | .671 379 | —.09426 94456 39 | .621 
.330 .45296 76932 02 .670 .380 | —.26632 63801 78 | .620 
.331 .41365 95766 99 | .669 .381 | —.48045 19086 93 | .619 
332 .52433 35362 36 | .668 | .382 | —.55645 23195 61| .618 
.333 .45999 71920 26 | .667 .383 | —.52637 52089 99| .617 
334 .37014 12343 31| .666 | .384 | —.33212 51324 29| .616 
335 .22552 02235 19 | .665 385 | —.39479 74417 87| .615 
.336 10904 55459 52 | .664 .386 | —.54657 07923 33| .614 
337. |  .23421 37763 08 | .663 .387 | —.66192 84709 16 | .613 
.338 .25303 23429 37| .662 | .388 | —.69040 53292 05 | .612 
.339 .23376 32937 51| .661 .389 | —.49978 37896 32| .611 
.340 .05089 90862 54 | .660 390 | —.48100 38625 93 .610 
341 | —.15716 90669 72 | .659 391 | —.50444 79258 90 | .609 
342 | —.09044 72949 99| .658 | .392 | —.56183 05832 68 | .608 
.343 | —.01752 94215 89 | .657 .393 | —.62318 58803 80 | .607 
344 .09699 75179 51 | .656 394 | —.50778 26196 53 | .606 
345 | —.01688 93672 23| .655 .395 | —.49647 33948 32 | .605 
346 | —.25654 09168 79 | .654 396 | —.41302 52891 40| .604 
.347 | —.27382 88597 64| .653 | ..397 | —.35589 92135 83 | .603 
.348 | —.21463 82850 53 | .652 | .398 | —.38867 82210 92| .602 
349 | .00184 25835 82| .651 .399 | —.35791 63481 26 | .601 
| —W(z) s -W(2) = 
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On the Numerical Solution of Convolution 
Integral Equations and Systems 
of such Equations 


By J. G. Jones 


1. Introduction. This paper discusses the application of a simple quadrature 
formula to the numerical solution of convolution integral equations of Volterra 
type and to systems of simultaneous equations of the same type. The convergence 
of the processes is considered in some detail, proofs being given that at a fixed 
value of the independent variable the errors in the solution tend to zero as the 
step length tends to zero. 


2. Types of Equations Considered. The convolution integral equation to be 
solved is 


(1) agit) = [ W(x — £)g(t) dé = f(z) 


where ‘a’ is a constant, and f, W are given functions. The cases a = 0 (equation 
of the first kind) and a ~ 0 (equation of the second kind) are discussed. 
The corresponding system of equations can be written 


(2) Ag(z) — [ W(x — &)g(t) dé = f(z) 


where 


Gn, °°*. Oe 
A=]: : |,a@ matrix of constants, 
Ani sites Onn 


Wulz)--- Win(z) filz) gi (x) 
War(z) eee Wan(z) f(x) 9n(x) 


and the notation involving integration is interpreted in an obvious way. For this 
system of equations the rank of A is an important parameter. 

Equation (1) is well known, commonly occurring in practical problems. The 
system of equations (2) does not appear to have been previously discussed, so an 
instance of a practical problem in which it arises is briefly described. 

In the linearized supersonic theory of [1] the system of equations (2) gives the 
pressure coefficients f;(z) on n spanwise wing stations produced by a quasi-cylindri- 
cal shaped fuselage defined by n Fourier components g;(1). The W;;(z) are tabu- 
lated influence functions and A is of rank ‘unity. In [2] the problem of determining 
the fuselage shape required to produce prescribed pressure coefficients on wing 
stations is considered. Special attention is there devoted to the numerical solution 
of (2) when n equals 2 and the rank of A is unity. 
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3. The Use of the Trapezoidal Rule. In the method to be described, the equa- 
tions are solved by dividing the z-axis into equal small intervals, replacing the 
integrals in the equations by the corresponding approximate expressions given 
by the trapezoidal rule, and hence deriving a formula for the step-by-step solution. 

The solution of an integral equation by replacing the integral involved by a 
quadrature formula may be regarded as a standard method; see, for example, [3]. 
However, adequate discussions of the convergence of such numerical methods 
appear to be lacking. This paper is concerned with the convergence of the method 
in the particular case of the application of the trapezoidal rule to equation (1), and 
with the extension of the method to the solution of the system of equations (2). 

The trapezoidal rule is the least accurate of quadrature formulas for a given 
step size, but it has the computational advantage of leading to a simple routine 
suitable for either desk or automatic computing, and furthermore is the most simple 
case to consider as regards convergence. 

Subject to assumed bounds on the derivatives of the functions involved, it has 
been found possible to establish a bound for the truncation error at a given value 
of the independent variable, and hence to prove the convergence of the process. 

The method of establishing a bound is not intended as a practical means of 
estimating errors. Practical ways of checking the accuracy of a numerical compu- 
tation of this type are either to repeat the solution using a different step length, 
or to evaluate the integrals in the equations by using the numerical solution and a 
more accurate quadrature formula than the trapezoidal rule. However, it is shown 
in Section 8 that for equations of the first kind the numerical method described may 
in some cases also provide a convenient means of error analysis. 


4. Sources of Error in the Solution. Consider the truncation error in this type 
of numerical solution. In going from the mth step to the (m + 1)th, errors occur 
from two sources (neglecting rounding errors which it is assumed throughout are 
kept smaller than the truncation errors) : 

(i) due to replacing the integrals by a quadrature formula 

(ii) due to the fact that the values of the solution at the first m steps, all of 
which go into the quadrature formula for the evaluation of the solution at the 
(m + 1)th step, have corresponding errors in them. 

In the following treatment the error at a fixed value of x is considered and the 
step length is taken as an integral fraction of x. As the step length decreases, the 
decrease in the corresponding error produced in a single step is offset by the fact 
that the number of steps required to reach x increases inversely. An analogous 
error analysis for the numerical solution of an ordinary differential equation is 
given in [4]. The form of the bound obtained in this latter case is similar to that of 
the bound described in Section 6. 


5. The Single Integral Equation. In this section we consider the numerical 
solution of (1). If a + 0 the step-by-step process is started by means of the equation 


(3) ag(0) = f(0). 


If a = 0, (1) gives, on differentiation, 
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(4) ~W(0)9(z) — [ W(x — &)g(&) dt = f"(z). 


If now W(0) ¥ 0, the step-by-step process is started by means of the equation 
(5) —W(0)g(0) = f’(0). 

If W(0) = 0, further differentiation (mn + 1) times is required until we reach a 
derivative W‘” (0) # 0. 


Equal intervals of length 6 are now taken along the z-axis. For the mth step we 
have 


mé 
ag(ms) — l W (ms — £)g(t) dt = f(ma). 


Replacing the integral by the corresponding trapezoidal rule approximation and 
rearranging : 


\a -3 w(0)}9*(ms) = f (ma) 
(6) 1 and 
+ if} W(ms)g(0) + > W({m — aagr(ia) } 


where an asterisk denotes an approximate value (and we use the convention 
throughout that > >m = 0). 
Unless both a = 0 and W(0) = 0 we can successively evaluate 


g*(5), g* (26), a 


If a = Oand W(0) = 0, then the equivalent equation (4) (or if necessary an equa- 
tion obtained by further differentiation) is solved by the method described. The 
additional errors involved in numerical differentiation (assuming the functions to 
have been given in tabulated form) are discussed in Section 8. 

The equation corresponding to (6) involving the true values is 


{a —§ w(0)} gms) = (ms) 


(7) _ ™ 
+6 ‘3 W(méb)g(0) + > W([m — asdo(ia)} + p> Cm, i 
where 
em = | Wb — e)g(t) at 
(i—1)6 

(8) ‘ 

— 55 Wi — [¢ — 1]6)g((i — 1)6) + Wms — is)9(ia)} " 
Writing 


(9) En = g*(mb) — g(ms) 
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(i.e., Z,, is the error after m steps) equations (6) and (7) give 
m—1 ™ 

(10) {<a - wo) En = 6D) W([m — 6) E; — Do ems 
i=1 t=1 


Now (see, for example, [5]) from (8) 


wo he 
12 


Cn,i = 


2{d d . \ 5 
aad — eee _ , 

6 (4 [W(mé — £)g(é)l mus di [W(mé £)9(E)eacsnay + O(8). 

If we assume that g and W have bounded derivatives up to the second, it follows 

that en: is O(8). 


From equation (10) it follows by induction that a positive constant K can be 
found such that 


| En | = fm 


where 
3 \ m—1 . 
(11) a — 5 WO) fm = OX fe + mK. 
t=] 
Rewriting (11) with (m + 1) for m and subtracting (11) from the result, 
\a _ ; W(0)} fan = {a — : W(0) + aK Sn = K®. 


The solution of this difference equation with initial condition from (11) is 





a— 2 w(0) + 0K | 
(12) fa = Ci | pevegi one Md 
{ ” 5 W(0) 


where 
Ké + (a a w(0)) 
2 
C; SS i 7 Soma inion a 
ot W(0) + 6K 


The two cases a ~ 0 and a = 0 are discussed separately in the next two sections. 


6. Equation of the Second Kind (a ¥ 0). In this case (12) provides the required 
bound. It is similar to the bound obtained for the error in the numerical solution of 
an ordinary differential equation in [4]. 

At a fixed value of x we consider a sequence of values of 6 such that mé = x and 
m— © as §—0. Then for small 6 the term in (12) with index m is bounded and 
C, = 0(8’). So for a fixed value of x, fm = 0(8'), hence E,, = O(8). 

By considering the particular case W = constant = K > O and g” = constant = 
12 it can be seen that the error E,, can attain its bound f,, , and so, for given step 
size, it is possible for E,, to increase with m like A”, where A > 1. 
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7. Equation of the First Kind (a = 0). We assume in this section that W(0) < 0, 
in which case equation (1), with a = 0, can be solved directly by means of the 
numerical formula (6). 


In this case equation (12) only provides a convergent bound for E,, (at a fixed 
value of zx) if 


(13) 


2 Se 


~ WO@|*” 


In general K cannot be chosen to satisfy this condition. An example of an equation 
for which (12) does provide a bound is 





(14) [ g(t) dt = 22°, 


in which W = constant = 1. In this case it is easily shown that the errors satisfy 
the equation 


m—1 
(15) — 5 En =D Est mi. 
On comparison with (11) it can be seen that we can choose K = 1 in order to satisfy 
the equality E,, = f,, . In this case (13) is satisfied, and as in Section 6 it follows 
that for a fixed value of z, E,, = O (&). A point of interest is that the errors E,, 
in this example form an oscillating sequence, viz., —28°, 0, —28°, 0, ---. This 
behavior is typical of the errors when a = 0, as is shown by the following analysis, 
which also shows that a convergent bound can always be found. 

Equation (10) in the present case reduces to 


(16) — 2 W(0) En = >> W([ m — il8)E, — 7 — 
So 

(17) ° 5 W(0) (Ems: — E,.} = 8W(8)En + Ju 
where 


m—1 
Jm = 85>, {W([m + 1 — i]s) — W([m — a8) }E; 
t=1 


(18) » 
— Cm+im+1 — > (€m41,i — Ons}. 


Equation (17) can be rearranged to give 


(19) = 


wlio 


® W(0) {Ems + En} = 8{W(5) — W(0)} Em + Jn - 

The right hand side of (19) is in general smaller than the right hand side of (17). 
That is, Eni; + E,» is in general smaller than E,,,,; — E,,, which implies that 
En4, and E,, are in general of opposite sign. The £,, then forms an oscillating 
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sequence. A smaller bound than that given by equation (12), which is obtained 
essentially from (10) by replacing the errors by their moduli, can therefore be 
obtained by relating E,,4. to E,,. 

Replacing m by m + 1 in (19) and subtracting (19) as it stands gives an equation 
for Ens2 — E,,. Assuming that g and W have bounded derivatives up to fourth 
order, it follows that the first and second differences of W are respectively of orders 
§ and & and that the first and second differences of ém,; are respectively of orders 
5‘ and 8. It can then be shown inductively that a positive constant K can be 
found such that 


(20) | Em | S hm 

where 
m—l1 

(21) Ima — Tim = Kis Do he + 8(Rmgt + hm) + mst + i, 
t=] 


and h, , he are chosen so that (20) holds for m = 1, 2. 
From (16) it can be seen that h; and h. can be chosen so that 


(22) kh, = O(8), he = O(8’) 
and 

(23) hi < he. 

Then 

(24) hom—1 < hem, ail _m 


and hence it may be verified by induction that if 


m—1 


(25) Hemsz — Hom = K {a(2 > Ae + Hm) + i( Han aa Hansa) + 2ms* + ‘ 


t=1 


and 


(26) Hz = he 
then 
(27) hem < Hom all m> 1. 


From (25) there is obtained the difference equation 
(28) Hemss(l — Ké) — Homae (2 + K&’) + Hom(1 — K&* + Kb) = 2K8*. 
The solution of (28) is 


(29) Hom = -s + Cyp." + C2p2” 
where p; , p2 are the two values of 
1 Ke’ 2 72.63 K*s* }\” 


and C, , C2 are given by the initial conditions 


sr + pC: = H, +8 


(31) 2 2 2 
Ci + pC, = Hy + 6. 
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Now from (22) and (26), Hz = O(8) and from (25), Hs = HA + O(8)). 
Also from (30), p, = 1 + O(8), pp = 1 + O(8), m2 — ~: = O(8). So from (31), 


c, = PH + ©) — (Ha + &) _ 068") 


ae = gens Bp 





Similarly C. = 0(8). 

At a fixed value. of z we consider a sequence of values of 5 such that 2mé = z, 
thus m — © as §— 0. Then for small 6 the terms in (29) with index m are bounded 
and C; , C2 are O(8'). So for fixed z, Hom is O(8) and from (27), (24) and (20) it 
follows that E, = O(8&). 


8. Comparison of Methods. An equation of the first kind (a = 0) can be con- 
verted into an equation of the second kind by differentiation. For example, if 
W(0) + 0 equation (4) results. It has already been shown in Sections 6 and 7 that 
for both types of equation the truncation error at a given value of the independent 
variable is 0(8°). In this section the effects of rounding errors in the two cases are 
briefly considered and the problem of whether to convert a first order equation 
into a second order one by differentiation before solving numerically is discussed. 

In the case a = 0 the numerical formula (6) becomes 


(32) gts) = pigs {2 — [5 wremayo(o) + Welm — aadgr(aa) 


In the first place, as was shown in Section 7, the truncation error in g*(mé) obtained 
by using the formula (32) is in general of opposite sign for two consecutive values 
of m. Because of this the truncation error in the summation in the right hand side 
of (32) and hence in g*(mé) is smaller than it would be if the errors were of the 
same sign. To take advantage of this in the numerical computation, several digits 
must be retained beyond the point where the truncation error begins before round- 
ing off. This implies that W must be known to several more significant figures than 
would otherwise be necessary. It is also evident on comparison of (6) and (32) 
that in the latter formula more digits in f(mé) must be retained before rounding 
because of the 6 in the denominator. 

Suppose now that we are to solve an equation of the first kind, viz. (1) with 
a = 0. The functions f and W are supposed to have been given in tabular form and 
W(0) =~ 0. Then the equation can be solved directly or it can be differentiated 
first, giving (4), an equation of the second kind. 

If the first method is adopted, then, as has been shown, the truncation errors 
will in general be of opposite sign and the rounding errors must be kept several 
digits smaller than the truncation error. If an automatic computer is being used 
it may not be inconvenient to retain these extra significant figures. However, f and 
W must be known to the extra degree of accuracy. Since the truncation errors are 
in general of opposite sign a simple smoothing process may be employed to improve 
the solution as follows. Denote the sequence of numbers 


g(0), g*(25), g*(46), --- 





by the symbols 
ge(0), ge( 25), ge( 45), --- 
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and complete the sequence 


ge(O), ge(8), ge( 26), + += 
by interpolation. 
Similarly, from the sequence of numbers 


g* (6), g* (36), g* (58), sh 
the sequence 


go(5), go(26), go(35), --- 
ig formed by using interpolation. 

The two sequences g.(mé4), go(mé) then give approximate bounds to the solution, 
and the smoothed solution is given by g,(mé) = 3[g.(mé) + go(mé)]. This procedure 
is illustrated in the numerical computation at the end of this section (Table 1). 

If the second method is adopted, that is, equation (4) is solved instead, the 
functions f and W must first be differentiated numerically. Once this has been done 
the numerical solution can be rounded off to the same degree of accuracy as the 
truncation error and consequently fewer significant figures need be retained. How- 
ever, extra significant figures in f and W now have to be used in the first place to 
obtain f’ and W’ to the required degree of accuracy when using numerical differen- 
tiation. 

Each method appears to have its advantages, and the choice must depend on the 
data provided in a given problem, that is, on the spacing and number of significant 
figures in the given functions. 

We conclude this section by presenting the details of a simple numerical computa- 
tion in which the results given by the two methods can be compared. 

The equation is 


[ “Wa — &)g(t) dt = f(z) 


where W(x) = cos x and f(x) = sin x. The analytical solution is g(¢) = 1. The 
functions are given to 4 significant figures at intervals of 0.1 in x. The functions are 
differentiated numerically using a three-point formula except for the values at 
x = 0 where a four-point formula is used. The interpolation in the formation of the 
sequences g.(mé), go(mé) is linear. 
































TABLE 1 
| | W'(x) | 
Ww) He) | Girect | gets) | gotey | sete | Groint | foe | cainbnat 

* | Goose) | (= sins) | Shuto | “| f° | | numerical | or Rula) asain) 

| 

| 
0 1.0000 0 1.000 sanry 1.000 | 1.000} 0 | 1.000 1.000 
0.1 0.9950 0.0999 1.003 | 1.000 | 1.003 | 1.002 | —0.100 | 0.994 0.999 
0.2 0.9800 0.1988 1.000 | 1.000 | 1.000 | 1.000 | —0.198 | 0.978 0.998 
0.3 0.9554 0.2954 0.997 1.003 | 0.997 | 1.900 —0.294 | 0.953 0.997 
0.4 0.9211 0.3894 1.006 | 1.006 | 0.998 | 1“902 | —0.389 | 0.920 | 0.998 
0.5 0.8776 0.4795 0.998 | 1.004 | 0.998 | 1.001 | —0.479 | 0.876 0.998 
0.6 0.8253 0.5647 1.003 1.003 | 0.996 | 1.000 | —0.564 | 0.823 | 0.997 
0.7 0.7649 0.6441 0.995 | 1.006 | 0.995 | 1.000 | —0.642 | 0.763 | 0.997 
0.8 0.6968 0.7173 | 1.008 | 1.008 | 0.994 | 1.001 | —0.716 | 0.696 | 0.998 
0.9 0.6216 0.7833 | 0.994 | 1.008 | 0.994 | 1.001 | —0.783 | 0.621 0.998 
1.0 | 0.5402 0.8415 1.009 | 1.009 | 
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9. The System of Integral Equations. In this section we consider the numerical 
solution of the system (2). We use the expression r(A) to denote the rank of A. 
If r(A) = n, the step-by-step process is started by solving the algebraic set of 
equations 


(33) Ag(0) = f(0). 


If r(A) < n then n linear combinations of the equations (2) can be chosen to 
give another system of the same form, with a matrix in which n — r rows are identi- 
cally zero. If the n — r corresponding equations are differentiated, the rank of the 
matrix of the resulting set of n equations will in general have increased (in the same 
way that differentiation of the equation of the first kind, viz. (1) with a = 0, in 
general gives equation (4) of the second kind). We assume that repeated applica- 
tion of this process eventually leads to a system of the same form as (2) with a 
matrix of rank n (this corresponds to the assumption in Section 5: W‘” (0) ~ 0 
for some n). Setting x = 0 in this new system, we obtain a set of equations anal- 
ogous to (33) with which to start the step-by-step process. 

Equal intervals of length 5 are now taken along the z-axis. For the mth step we 
have 


mb 
Ag(mé) -{ W(ms — £)g(E) dé = £(ma). 


Replacing the integral by the corresponding trapezoidal rule approximation and 
rearranging: 


{A +d Swio)} g*(ms) = £(ms) 
(34) m—1 
+ 645 Wima)g(0) +S Wilm — sa)grtia)} 


= Q(ms) (say), 
where an asterisk denotes an approximate value. 


Provided | A — * w(0) | # O we can successively evaluate 


g*(5), g*(26), --- 


where at each step we solve a set of n equations (with coefficients independent of n). 

The above method is exactly analogous to that described for the single integral 
equation in Section 5. Only the outlines of proofs of convergence are given in this 
section, the details being analogous to those given for the single integral equation. 
The equation for the truncation error corresponding to (10) is 


(35) {A ~ * wio)} E,, = >> ‘W({m — iJ6)E, — > en. 


where 


E,, = g*(md) — g(md) 
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and 
16 
enc = | Wms — s)g(e) dt 
(i186 
— 5 (Wms — [i — 18)g(fi — 118) + W(md — i8)g(68)}. 


The case of r(A) = n corresponds to the case a ~ 0 in Section 5 and a quantity 
fm can be found which is a bound for every component of E,, and is of the same form 
as (12) (with a ~ 0). So when r(A) = n, at a fixed value of x we have E,, = 
O(é°). In this case there is no condition on r(W(0) ). 

If r(A) < n the analysis proceeds as in Section 7. In this case we assume, to 
avoid further complication, r(W(0)) = n. The equation corresponding to (17) is 


(36) {A “ Swio)} {Enis — En} = 8W(8)Em + Jn 
where 

Jn = 8 > {W([m +1 — a8) — W(lm — 18) }E, 
(37) t=1 


m 
— Cntim+ — >» (@m4-1, — @n,:). 
i=l 


Equation (36) can be rearranged to give 
(38) {A - £wio)} {Eni + En} = 2AE,, + 5{W(5) — W(O)}E, + Jn. 


This is of the same form as equation (19) except for the term 2AE,,. It is now 
shown that when r(A) < n this term only contributes terms to the elements of 
E..4: + E, of the same order as the contribution oi the remaining terms on the 
right hand side of (38). We assume that the elements of the last n — r rows of A 
are all zero (this can be arranged by taking linear combinations of the original 
equations). Then 


rye wo) | = 0(8"”), 


and it is easily verified that each element of E,.,; + E,. consists of a linear combina- 
tion of the elements of the right hand side of (38), the coefficients of the first r 
elements being 0(1) and those of the last n — r elements being 0(3~*). But all the 
last n — relements of AE,, vanish, so the dominant terms in each element of En4; + 
E,, are of the type 


(a) O(8*){a[W(5) — W(0)]En + Im}: 
where r+1Sisn, 
(b) O(1) {AE,,}; 


where 1 Si Sr. 


The term (a) is analogous to the expression for En», + E,, from equation (19). 
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From equation (35) 
m—1 ™ 
"ies 5 W(0)En + 8D Wm — i8)E, — 3 en.s. 


It can be seen that contributions from (a) and (b) are of the same order. 

The situation is the same as in Section 7, the elements of E,,,, + E,. being in 
general smaller than those of E,.,; — E,, , consecutive errors thus being in general 
of opposite sign. Replacing m by m + 1 in (38) and subtracting (38) as it stands 
gives an equation for E,,,2. — E,,. Assuming that the elements of g and w have 
bounded derivatives up to the fourth order it can be shown inductively that a 
positive constant K can be found such that 


\{E,j:| She, foralli 


where h,, is given by (21), and (22), (23) hold. The analysis of Section 7 then 
shows that for a fixed value of x we have E,, = O(8). 

As indicated at the beginning of this section, if r(A) < n the system can in 
general be converted into an equivalent system with r = n by application of the 
operations of addition and differentiation. The problem of whether to convert a 
system with r < n into an equivalent system in this way, before solving numeri- 
cally, is exactly analogous to that discussed in Section 8, and will not be treated 
further here. The merits of the two methods are summarized in the conclusions. 


10. Conclusions. The numerical solution of a convolution integral equation of the 
Volterra type, equation (1), by using a simple quadrature formula has been dis- 
cussed. It is shown that if the step length is 6, the truncation error at a fixed value 
of the independent variable is 0(6’) both for equations of the first and second kinds. 
However, the behavior of the truncation error is different in the two cases, being in 
general of opposite sign for consecutive steps in the case of an equation of the first 
kind. Since an equation of the first kind can in general be converted into one of the 
second kind by differentiation, it can either be solved numerically as it stands or 
differentiated first. The merits of the two methods appear to be as follows. In the 
direct method, both a smooth solution and approximate bounds for the truncation 
error can be obtained simultaneously. However, more significant figures have to be 
retained throughout before rounding off. If the equation is differentiated first, the 
given functions (which are assumed. to be given in tabulated form) have to be 
differentiated numerically. In the actual step-by-step computation fewer significant 
figures need be retained. If the given functions are tabulated to sufficient significant 
figures to permit accurate derivatives to be obtained using numerical differentiation, 
the second method will probably be the more accurate. Otherwise, and if it is 
convenient to retain the extra significant figures throughout, the direct method 
has the advantages of requiring fewer steps (no numerical differentiation) and of 
providing approximate bounds for the truncation error. The method has been 
generalized to include the numerical solution of the system of integral equations 
(2). An instance of a practical problem arising from supersonic linearized theory in 
which this equation occurs has been described. The truncation error at a fixed value 
of the independent variable is again 0(8’). The behavior of the truncation error is 
analogous to that in the solution of the single integral equation, the case r(A) < n 
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corresponding to an equation of the first kind, and the case r(A) = n corresponding 
to an equation of the second kind. 
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Numerical Integration Formulas of Degree 3 
for Product Regions and Cones 


By A. H. Stroud 


1. Introduction. Here we discuss approximate integration formulas of degree 3; 
that is, formulas which are exact for polynomials of degree < 3, of the form 


[ sa, +++, ay) dt, +++ dti™ D> bi f(va, --+, vin) 


for certain regions R in n-dimensional, real, Euclidean space £,, . 

If R is the Cartesian product of an r-dimensional and an s-dimensional region, 
R = R, X R,,7r + 8s = n, and if formulas of degree k involving p and q points are 
known for R, and R, , then Hammer and Wymore [4] have shown that a formula of 
degree k may be obtained for R which involves pq points. Similarly, if R is a cone 
(or pyramid) with (n — 1)-dimensional base B, and if a formula of degree k involv- 
ing p points is known for B, and further, if a formula of degree k of the form 


[ a” g(x) dx ~ > aig(v:) 


is known involving qg points, then Hammer, Marlowe and Stroud [2] have shown 
that a formula of degree k can be obtained for R which involves pq points. 

For k = 3 we will show that these results can be improved as follows. If R = 
R. X R, and formulas of degree 3 involving p and gq points are known for R, and R,, 
then a formula of degree 3 can be found for R involving p + q + 1 points (and in 
some cases this may be reduced to p + q or p + q — 1 points). If Ris a cone and a 
formula of degree 3 involving p points is known for its base B, then a formula of 
degree 3 involving p + 3 points (in some cases p + 2 or p + 1 points) can be 
found for R, whereas the method of [2] would involve 2p points. We also give a 
p + 2 point formula of degree 3 for certain double cones, where p is the number of 
points in a formula for the base, and a 2n + 3 point formula of degree 3 for an 
n-dimensional simplex. 


2. Formulas for Product Regions. For convenience let us assume that the 
centroid of R, is at the origin of coordinates in the subspace E, of E, for which 
Xr41 = +++ = XL, = 0, and also that the centroid of 2, is at the origin in the subspace 
E, for which x; = --- = x, = 0. Then the centroid of R is at the origin in £, . 
Also, for convenience in notation, we will write 


Ray" +++ 2,7) = / ai +++ a,” dz, +++ ax 
R, 
with similar notations for R, and R. Then, for example, R,(1) is the r-dimensional 
content or volume of R, and R(1) = R,(1)R,(1) is the n-dimensional volume of R. 
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Suppose an approximate integration formula of degree 3 for R, is given by: 
(1) Vi = (Va, °°, Vr a;, i=1,2,---,p. 
Also suppose that a formula of degree 3 for R, is given by: 
(2) pti = (Mptirtis *** > Mptin) Ap+j, j = 1,2,---,¢. 


At first it is assumed no y; coincides with the origin in EL, and no vy; coincides with 
the origin in E, . An integration formula for degree 3 for R is given by: 


vw = (0,---,0,0,--- , 0) bo = —R(1) 
(3) vi = (va,°°+ , Vir, O,--- , 0) b; = a.R,(1), 
To Seren 
pti = (0, --- , 0, Mptirtts °° * » Yrtin) bp+j5 = Gp+jR,(1), 
j=1,2,---,¢@. 


To prove this statement, it suffices to prove that the formula is exact for the follow- 
ing ten types of monomials: 1, Z,, , Xs, 5 LryLrq 5 LryLe, » La, Leg » LryLryLry y LryLryLe, y Lr, Ley 
Log» Ley LeqXo, Where 7; , 72,73 = 1,2, --- , rand s, 8,3 =r+i,r+2,---,n: 

The monomial 1. Since a; + --- + a, = R,(1) and ay4; + --- + apy, = R,(1) 
we have 
—R(1) + far + +++ + ap|R.(1) + [apy +--+ + Ap+q)R-(1) 

= —R(1) + 2R,(1)R,(1) 


The type 2,,%,,7,, . The formula gives for the integral of z,,z,,x,, over R 


R(1). 


II 


[Qi 1r,V1rg¥ir; + e104 + OV pr\V pryV pra) s(1) = R(X r,%r.Xr,) Rel 1). 
Thus the formula is exact because 
R(X, %r—hry) = Rr(Ar,%r_Lr,)R,(1). 


The proofs for types 2,, , Xs, , UryVrq 5 Le,Xeq » Lo,Xe_%s, are similar. 
The type 2;,2,,2., . By the formula the integral of z,,z,,x,, over R is zero since each 
term of the sum is zero. However 


R(27,Xr_%e,) = Re(Xr,%r_)Rs(Xs,) = O 


because R,(z,,) = 0. The proofs for the types z,,x,, and 2,,%,,2;, are similar. 

If the formulas for R, and R, already include the origin, then the formula for R 
will involve either p + g or p + gq — 1 points according as the origin is included in 
one or both of these formulas. In the former case, if vp; in (2) is the origin, then 
bo = —[Gp42 + --- + Gpy_]R-(1). In the latter case, if both » in (1) and v4, in 
(2) are the origin, then 


bo = R(1) — [ae + +++ + ap|R(1) — [apse + ++ + GpyclR-(1). 


Depending on the particular structure of the formulas (1) and (2) it may be 
possible in special cases to eliminate the origin in the formula (3) for R. Suppose 
in the formula (2) for R, , vp4: is the origin and a,,; is positive. Also suppose R, 
is centrally symmetric. Then, as shown in [6], we may obtain 2r point formulas for 
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Ry with — 745 = vi, Gr4i = G;,% = 1, 2,---, 7. In addition to the assumptions 


already made about R, , we may also assumeR,(z,) = --- = R,(z,7) andRAza;) = 
0, for i + 7; then in order that (1) be a formula of degree 3 it is necessary and suffi- 
2 
cient that »,---, v, be orthogonal vectors with a; = ee where | »; | is the 
vi 


distance of v; from the origin. Since in the formula for R, X R, it is no longer neces- 
sary that b; + --- + be = R(1), a formula of degree 3 for R, X R, is given by 
the 2r + q — 1 points: 





—r45 = 15 = (a, --- , vr, 0, --- , 0) bres = 1, 
t= 1,2,---,r 
Vorsj = (0, --- , 0, vorsirtt, “°° » Martin) bers 5 = Ger+;R(1), 
j=2,---,@ 
where the only conditions that the b; and »;, i = 1, 2, --- , r, must satisfy are 
2(b: + -+- +5.) = GeryiR-(1) and »,---, », are any set of orthogonal vectors 
for which | »;|? = = ") ; the b; must all be positive. 


Of course if both R, and R, are centrally symmetric then R, X R, is also centrally 
symmetric and the results of [6] may be applied directly to obtain 2n point formulas. 
Some of the resulting formulas may also be obtained from the separate formulas 
for R, and R, by a method somewhat similar to that described in the preceding 
paragraph. 

As an example of a specific formula we give a formula for the region C, X S,_,, 


1 <=rasn -— 2, where C, is the r-cube with vertices (+1, +1, --- , +1) and S,_, 
is the (n — r)-simplex with vertices (—1, —1,---, —1), (1, 0,---, 90), 
(0, 1,---,0),---, (0,0,---, 1). For r = 1, C, X S,. can be considered as a 


prism (or in the terminology of Sommerville [5], p. 113, a prism of the first species) 
with base S,_; ; for general r, C, X S,_, is a prism of species r. Since C, is centrally 
symmetric the remarks made above apply. A formula of degree 3 is known for 
S,-, using n — r + 2 points, [3], one of which is the origin, but the origin is taken 
with a negative weight so that it is necessary to use n + r + 2 points in the formula 
for C, X S,.,. A particular formula is: 


vw = (0,0, ---,0,0,0, ---,0) = 


-7-Gorte) 
= 800 [1 5 Sear em 
—v41 = 1 = (1,0, ---, 0,0, 0, ---, 0) bar = bh = +--+ = 


r—1 


— Vr+2 == (0,1, --- ,0,00,--- , 0) = bh, = = 28) 
—v, = v, = (0,0, ---, 1, 0,0, ---, 0) 
—2 
Yr = (0,0, "nae =—r+3’ 
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—2 —2 
quamitneadieeeeee ee _ => ee = t = 
n—r+ 3’ "n—rt :) Darts Peter 


ye 2 _ (n—1r+3)? 28,_(1) 
mere (O05 GF go 0) = geet se 





2 
Vn+r+1 = (0, 0, ee, 0, 0, 0, ‘Ag — 5) 


3. Formulas for Cones. Suppose the centroid of the (n — 1)-dimensional base B 
is at the origin in the subspace F,_, of E, for which 2; = 0 and that the vertex of R 
is at (1, 0, --- , 0) in Z, . Suppose further a formula of degree 3 for B is given by 


(4) Vi = (V2, Vis, *** 5 Vin) aj, = 1,2, ---,p. 


We assume at first that no »; coincides with the origin in E,_,. An integration 
formula of degree 3 for R can be found of the form 


om (n, YVi2, YHi8, * °° » Vin) b; = ad;, i= 1,2, °-- » P 
Yeij = (§;,0,0,--- , 0) b;, j = 1,2, 3. 


To prove this we must first calculate the monomial integrals over R in terms of 
those over B. These are 


(5) 


[ ayaa +++ ay” dat, +++ day 
R 
1 
= [ (1 — 2)" thet tFa1 8 dar, / ae”? +++ 2," dr, +++ dite. 
0 B 
The equations which must be satisfied by (5) if it is to integrate each monomial are 


given below (where the monomial from which the equation arises is indicated at 
the left): 


'(n — 1)! 
[x1"] b a + be bP + bs a8 + an’ a; + +++ +a,] = a B(1), 
6 = 6,1,2;3 
[x"xJ ayn lain: + -++ + apr] = pee, B(x;) 
(n+ 6+ 1)! ' 
B= 0,1,2 
Bi(n + 1)! 
[ay’x;2;] ayn’ lar; 1; ste ApVpi Vp) = (n+ B + )1 B(x;2;), 
6B = 0,1 
1 
[52 2x] aya VziVizVu +eee t+ Ap VpiV pj Vp) SS B(2;2x;2x) 
n+3 
where i, j, k = 2, 3, --- , n. Because we have assumed B(x;) = 0, the three equa- 


. 8 “ . te 
tions [z; z;] are of the form 0 = 0. From the equations [x,2;| and [x,7,x,;| it imme- 





oe fo os LI 


of 


re 
ut 





NUMERICAL INTEGRATION FORMULAS OF DEGREE 3 147 


diately results that 7 = 1/(n + 3). From equations [z,z,] and [z,7;) if follows that 
(n + 3)’ n+2 
~ (+2 an FB" 
Since a; + --- + a, = B(1) the equations that remain become 
[x,') at + baat + bet + OFS ay = OO pay, 
B = 0, 1, 2, 3. 


In these equations we cannot take, say, 6, = 0, which would in effect reduce by one 
the number of points in the formula, because the resulting equations do not have a 
real solution. It might be expected, however, that there would be many solutions 
with all b; + 0. Here we give just one of the simplest solutions. If we choose 


_+3) Eiy 
(n re ay BU), §3 _ n+ 3 





b; = 








then the other values are 


b/B(1) = 2(n + 1) — (n — 1) V2(n + 1)( n+ 2) 











b= 4n(n + 1? B(1) 
g, = 2m +2) + V2(n + 1)(n + 2) 
; (n + 2)(n + 3) 

be = by'B(1) = 2M +L) + (n — DV2An + 1)(m + 2) By) 





4n(n + 1)? 

2(n + 2) —+/2(n + 1)(n + 2) 
(n + 2)(n + 3) 
These values of b;’, be’, & , & have been given in connection with numerical integra- 
tion with respect to a weight function «"’ by Fishman [1] for n = 1(1)6 to 12D 
and in [2] for n = 2(1)4 to 18S. These authors have used n where we use n — 1 
and for m = 2,7 = 1, 2, their 1 — z; is our £;, their 6; is our b/. It is proved in 
[2] that the 2p points 
(&;, (1 — &))ve2, (1 — &))va,---, I £5) in) b,’a; 

a= 1,2,---,p,j = 1,2, are an integration formula of degree 3 for R. 

If one of the v; in (4), say » , is the origin in Z,_, then it may be possible to 
determine 6, , be, & , & with b; = 0 to give a formula which involves only p + 1 
points v2, --- , ¥p42. The formula using n + 2 points when R is an n-simplex, [3], 
can be derived in this manner; the formula for B, which is an (nm — 1)-simplex, 
involves n + 1 points of which one is the origin in Z,_, . In any event, if » is the 
origin, we may derive a formula of degree 3 for R involving p + 2 points », 
Vp+3 Where 


& = 





_(n+ 3)’; st 
~ (n + 238 


and the.other values as before. Usually there will be other p + 2 point formulas 
as well. 


+ +++ + ay) 
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Now we briefly discuss formulas for regions which are double cones; that is, 
regions which are the union of two cones with vertices (1, 0,---, 0) and 
(—1, 0, --- , 0) and with a common base B of the type we have just considered. 
By a combination of the methods of this section and the preceding it is easy to see 
that from a formula (4) for B we may obtain the following formula involving 
p + 2 points for the double cone R: 


ee Pyese — n+ 3)" sa 
i= (0, we, a; » Win) b; — “(n + 2)* ai, 7 > 1, 2, »Pp 
8 
—~Vore2 = Von = (é, 0, 0, -++,0) bois beri _ egy BU) 


where 











_n+2 se 2n? + 8n + 8 
1 ae8 one gm 4/ ett 8 


Since ¢ < 1 for all n, the points »,,; and v,2 are always interior to R. In this 
formula the origin is not required if it does not occur in (4) and b,4; and & are 
uniquely determined. If the origin does occur in (4), then in some cases it can be 
eliminated in the formula for R; if », is the origin, then 6,4; , &, and b, will not be 
uniquely determined and we may take ¢ arbitrary, which will then determine 


2 


bot = Bain F Im + D) POY 





2 
by = 2 BCL) — 2p, — M23) [a2 + +++ + ap). 


(n + 2)* 
Since 6,4: is required to be positive, we could, providing 
(n + 2)° 
[a2 + eee aa ay) <= nin + 3) B(1), 


or, in other words, providing 


_ (8n + 8) 


ay > n(n + 3) B(1), 
choose b, = 0. In this case the formula for R would involve only the p + 1 points 
vo, *** y Mpy2 and by; and — would again be uniquely determined. In this paragraph 


R has been a double cone of the first species with base B; this method may be 
repeated to give a formula for a double cone of species r. If we take a double cone 
of species r with base S,_,, 1 S r S n — 2, then this region has a formula of degree 
3 which involves n + r + 2 points. 


4. A Special Formula for Simplexes. We digress here from the methods of the 
previous sections to give a special formula of degree 3 for a simplex S,, in E, . The 
formula involves 2n + 3 points, n 2 2, all but one of which are on the surface of 
S, ; the formula involves the n + 1 vertices, the n + 1 centroids of the (n — 1)- 
dimensional faces, and the centroid of S, . 
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To develop the formula it is most convenient to take S, to have vertices 


(0, 0,--- , 0), (1, 0, ---, 0), ---, (0, 0,--- , 1). Then the monomial integrals 
over S,, are 





-la:! ! 

. . Ag+ Aj ay: 
[ ettapine® dey +++ dry = —— 
Sp (n + a; + aj + oy)! 


The formula is then: 


n= ( OB him =) by ow OU SS ie 





n+l1’n+1 "n+1 (n + 2)(n + 3) 
vy, = (0,0, ---,0) bh = +++ = bay = 
y= (1, 0, ---,0) : S,(1) 





~ (n+ D@ + Dn + 3) 


aH = (0, 0, ee, 1) 





ts eee Segoe oa 
1 1 in n° 
n43 = (0, = +4) ~ (n+ I(n.+ 2)m F+ 3) S.(1) 


3 
V2n+2 -(2.4, --+,0) 
nn 


where S,(1) = 1/(m!) is the volume of S,. Because of the symmetries of this 
particular simplex, the proof that this is a formula of degree 3 can be established 
by verifying that it is exact for the 7 monomials: 1, 2, , 2:°, x22, 21°, 21°22 » MXeks. 
For n = 1 this formula gives us Simpson’s formula, since for S, , a line segment, 
the vertices coincide with the (n — 1)-faces; the weight for each end point of 8, is 
b; + be = be + 63. For n = 3, bo = O, and in this one case the centroid does not 
occur in the formula; for n > 3, bo is negative. 

This formula will be useful when it is desired to integrate over a region by sub- 
dividing the region into simplexes and then applying a formula of degree 3 to each 
simplex. For a large number of subdivisions the total number of points used in 
applying this formula to each simplex will be less than the total number of points 
used if the n + 2 point formula is applied to each simplex. For example, for n = 3 
suppose we desire to integrate over an icosahedron (which has 20 triangular faces) 
by subdividing it into 20 tetrahedra (each tetrahedron having as vertices the 
vertices of a face plus the center of the icosahedron). Repeated use of the formula 
given here would involve a total of 63 points whereas use of the 5-point formula 
for each tetrahedron would involve 100 points. 


5. Concluding Remarks. From the formulas we have discussed, formulas may 
be obtained for regions which are linear transforms of the particular regions we have 
considered. While these results add greatly to our knowledge of formulas of degree 3 
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nothing is yet known concerning such formulas for regions which are less regular 
than those considered here. 

It seems certain, although no proof is known, that the n + 2 point formula of 
degree 3 for S, and the 2n point formulas of degree 3 for C, involve the minimal 
number of points for this degree. If this is true, then it is also likely that the n + 
r + 2 point formulas for both C, x S,_, and the double cone of species r with base 
S,_, are also minimal. For fixed r these latter regions are the duals of each other 
(see [5], p. 56) and thus have the same symmetries. This seems to indicate that the 
minimal point formulas of degree 3 for a region are related to the group of sym- 
metries of the region. 

I am indebted to R. G. Hetherington for many discussions on this subject. 
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The Epsilon Algorithm and Operational Formulas 
of Numerical Analysis 


By P. Wynn 


1. Introduction. It is the purpose of this paper to describe a non-linear technique 
which appears to have powerful and general application in numerical analysis. 
However, before doing so it is necessary to refer to a few related theoretical con- 
cepts. 


2. Rational Operational Formulas. The double sequence of rational functions 


mo u,v = 0,1,--- 
where 
(1) U,»(z) _ Ayo + Apri t + “oo a Ap 2” 





Verlt) — Bevo + Brat +--+ + Barn D 
may be derived from the series 
(2) B(x) = Dean" 


by imposing the condition that the power series expansion of (1) should agree 
with (2) as far as the term in z**’. If none of the Hankel determinants 


Cus Came °** Cm+k—-l 
Cm+1 Cm42°°° Cm+k 
realli wot Ge oe eek ee 
° 
Cn+k—1 Cmt+k °°°* Cm+2k—2 | 


vanish, and the additional condition 8,,,. = 1 is imposed, the coefficients in the 
rational expression (1) are uniquely determined. The rational expressions (1) may 
be placed in a two-dimensional array in which the quotient (1) occurs at the inter- 
section of the (u + 1)th row and the (v + 1)th column. [1] [2] [3]. 

As is well known, the numerical convergence of the sequence 


U,.(x) 
(3) V,2(x) 
for a particular value of x is in many cases much better than that of the series (2). 


This consideration led Kopal [4] to the consideration of rational operational forma- 
las, that is, to the replacement of the operational equation 


(4) (Sea)r=s 





r=0,1,--- 


s=(0 


where F is a known function from which f is to be determined, and d is a finite dis- 
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placement operator, by the equation 
vd)) 7 
ve Gan Tori 
Equation (5) cannot at the moment, in its non-linear form, be solved. The equation 
may however be linearized by multiplication throughout by V,,-(d) to give 


(6) U,,.(d)F = V,,-(d)f. 


Assuming that F and f are completely known. that r in equation (6) is sufficiently 
large, and the example is a suitable one, then there will exis‘ considerable numerical 
agreement between the right and left hand sides of equation (6). Assuming that d is 
any one of the conventional operators A, E, V, u, 6 of numerical analysis, and that 
F and the sequences of values fi , fo, --- ;f,f-2,-°-: are known, then equation 
(6) may be rearranged so as to determine fy . It is this very last assumption which 
constitutes a serious limitation of the linearizing technique resulting in equation 
(6). Indeed, Kopal was only able to find useful application of the technique when 
d was the backward difference operator, though his numerical results, which related 
to the forward integration of a differential equation, appeared to be very promising. 
However, the sameeffect over a very much larger range of problems may be achieved 
by recourse to another method. 





3. The e,,(S,) Transformation. In his researches into the non-linear trans- 
formation* 








Sn Sna4i a Snim 
AS,, ASn41 ria ASain 
: : oo : 
(7) €m( Sn) = ASn+m-1 ASn+m ii ASn+20-1 | mn = 0, 1, yee 
1 1 oe 1 | 
AS,, ASnsi 24> ASnim 
— ae. Seabee See 
of the sequence S,, r = 0, 1, --- Shanks [5], by an appeal to the theory of linear 
equations, showed that if 
(8) S,= > cx r=0,1,--- 
e=0 
then 
=n Un .min(X) wid 
(9) €m( Sn) Von.maa(Z) mn = 0, i‘. . 


The same result may be derived from the theory of orthogonal polynomials [6]. 


4. The e-Algorithm. The evaluation of the determinants in the various expres- 
sions (7) is sufficiently laborious to be prohibitive. However, the expressions (7) 


* The notation used here is consistent with that of [7] but differs slightly from that of 
[5] where the right hand side of (7) would be designated as ém(Snim). 
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may be computed recursively by means of the «Algorithm as follows I. If, from 
the initial conditions 


(102) &P=0 m=1,2,---; eo” =S, m=0,1,- 
quantities «s” are computed recursively using the relation 
m 1 
(11) en =e) + ape m,s = 0,1,---, 
€, _ 
then 
(12) Ath = fe(ASa))* fs = e(Sm) m8 = 0,1, - 
If the quantities «{” are arranged in the scheme 
Pts 
,: eee 
(1) 
€9 ° 
) ei” ” a 
a ole ey 
(3) (2) *. ight 
€1 €i €s 
(3) * (1) 
€0 F €s+1 
Pa 
(2) 


€s+1 


it will be seen that relations (11) may be used, column by column, to build up the 
scheme from left to right. It should be noted that if conformity, by means of equa- 
tions (9) and (12), is to take place between the Padé Table and the e-array, the 
latter must be transposed about the diagonal m = 0; the columns of the «array 
with even order suffixes then take their place as rows in the Padé Table. 

The following theorem, based upon the results of the last two sections, may now 
be given: 

THEOREM. /f p is an associative and commutative operator, and 





(13) a.p'F = c,x" s=0,1,--- 
and quantities es” are computed using the relation (11) from the initial values 
(14) & =0 m= 1,2,---; a” = > a.p'F 
s=0 
then 
- (m) __ Us mis(%) a4 vt 
(15) a ae Ve.mss(z) m, & 0, ' 


First ExamMpLe: A numerical example of the application of the theorem now 
follows. It concerns the process of obtaining the derivative at z = 0 of the function 
exp(hz), when h = 0.6, by means of the formula 


d _ = (—1)° pos 
(16) (gr) -5o% F. 








154 P. WYNN 


Here, in the notation of equation (13) 








(—1)° om . (78) +5 oti 
(17 ) c+4 Paap ¢ —) } 


The quantities with even order suffix in the ¢-array for this example are displayed 
in Table I. 

Note: The results of Table I begin with the diagonal m = 1. If the notation of 
equation (15) is strictly to be adhered to, an entry «)” = 0 together with a corre- 
sponding diagonal should be appended to Table I. However, this is not a matter of 
great importance, and in the event that an operational series were to begin with a 
term in p’,s > 1, even this artifice would not be available. 

It is perhaps in order to comment upon the power of the algorithm as revealed 
by this example. Attainment of the same accuracy as is achieved in Table I by 
the straightforward use of the series (16), even neglecting the accumulation of 
round-off errors, would involve the summation of about eighty terms and an excur- 
sion into arithmetic involving twenty-eight decimal figures. 

The Padé quotients (3) in this example are successive convergents of the con- 
tinued fraction 

= 1 Wx 1x Qa 2x 
(18) £ log (l+2) = —orarapat 
Numerical investigation into the behavior of this continued fraction [8] shows that 
application of the «algorithm to the series (16) converges quite reasonably for 
(e’* — 1) > 1, when the series rapidly diverges. 

In the derivation of the classical operational formulas of numerical analysis 
the operand is assumed to be a polynomial, and the formulas derived are then 
completely valid. The formulas are then universally applied, without examination 
of the operand, and without any more justification than that of the results achieved. 

In the same way it occurs that although formula (13) is no longer valid, use of 
the ealgorithm in conjunction with operational series meets with success. Two 
examples which support this thesis now follow. 

SeconpD EXAmpPte: This concerns the interpolation of the function log (0.6 + hz) 
when h = 0.1 and z = 0.25 with points of tabulation at unit intervals of z, by use 





0.8221 1880) 

0.4841 7914) 0.6038 2270) 

0.6693 9684) 0.5987 5229) 0.6000 7869 

0.5551 9362) 0.6005 0406) 0.5999 7937) 0.6000 0168 

0.6303 0451) 0.5997 6720; 0.6000 0665) 0.5999 9961) 0.6000 0004 

0.5788 4612) 0.6001 1785) 0.5999 9752) 0.6000 0011) 0.5999 9999}0.6000 0000 
0.6151 0747) 0.5999 3622) 0.6000 0102; 0.5999 9996; 0.6000 0000 
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0 
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of Bessel’s Interpolation Formula 


(19) F(z) = Sa F(0) 


= 


where 


a = 1 a, = 2k” 
_({2+8s-1 1/2 
a, = ( 2s ) aE , 


_ 8-§ 248-1 +1 ple, ae 
ett = oo 1 28 )e E's = 1,2, 


(20) 





The quantities «$”” with even order suffix are displayed in Table II. 


Since log(0.625) = —0.4700 036, it will be seen that application of the ~algorithm 
results in an effective gain of three decimal figures. This is not spectacular, but 
there is no point in selecting for presentation only those examples which display 
the method in a particularly favorable light. It might be mentioned at this point 
that the author has experimented with the e-algorithm in conjunction with opera- 
tional formulas in a large number of cases, and in none of these was the accuracy of 
the transformed results worse than the original partial sums. 

Since the odd and even order terms in the series (19) are so dissimilar, the odd 
and even terms were separated out and the two series submitted separately to 
treatment by the ec-algorithm, the transformed results subsequently being added 
together. The numerical results produced in this way were not, however, signifi- 
cantly better than those shown in Table II. 


TuirD EXxamMpPte: This concerns the application of the Euler-Maclaurin integra- 
tion formula 


m a? ™ 


> Gayl A,F™ (x) 





z+h h 
(21) [PW =a = 5 (F@) +F@ +h} - 


when the integrand is the function exp(—z’) and the upper and lower limits of 
integration are 0 and w(1 + 7) respectively, with w = 0.75. 
te! ; ;, 
The functions u, = ————~— in this example satisfy the recursion 


(s+ 1)! 
(22) s(s + 1)u, + 2h’su,. + 2h?(s — 1)u,2 = 0 


TABLE II 





—0.4337 503 
0.4722 880 | —0.4701 290 
0.4700 009 | 0.4699 404 | —0.4699 923 
0.4699 419 | 0.4699 732 | 0.4700 085 | —0.4700 040 
0.4700 084 | 0.4700 105 | 0.4700 032 | 0.4700 027 | —0.4700 034 








0.4700 105 | 0.4700 088 0.4700 027 | —0.4700 031 | 
0.4700 020 | 0.4700 017 0.4700 040 
0.4700 017 | —0.4700 019 


—0.4700 048 
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with 


—w{cos(2u*) + sin(2w*)} — iw{cos(2w*) — sin(2w*)} 


(23) 3 2 * 2 3 2 : 
—2w'{cos(2w*) — sin(2w*)} + i2w'{cos(2w*) + sin(2w*)}. 


Uy 

The quantities <{” (which are now complex numbers) with even order suffix 
are displayed in Table III. Since erf(0.75(1 + i)) = 0.917 306 + 10.403 654, 
application of the e-algorithm has in this case resulted in the gain of about three 
decimal places. 

It is perhaps of interest to point out that the accuracy of the transformed 
results produced in the first and third examples could have been increased by ex- 
tending the computation. This is also true to a limited extent of the second example, 
but the non-existence of central differences above a certain order limits the extent 
to which the computation may be prolonged. 

It would be useful, when examining the mathematical validity of the procedures 
adopted in the second and third examples, to be able to relate the determinantal 
quotient 





| s=0 s=1 s=m 
opr Lepr --- LeapF 
s=0 s=0 s=0 
c; pF Co pF “++ Oma pF 
| Ce p"F —_ p™ FP a Com p"F 
| 4 1 pen 1 
c; pF co pF -** Oma PF 
Con p"F Cots pF ee Com pF 
to the solution f of the operational equation 
p cp F - I 
s=0 


but this appears to be one of the cases in which a statement of the problem is not a 
great step forward to its solution. 
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Expansion of Spherical Bessel Functions in a 
Series of Chebyshev Polynomials 


By A. M. Arthars and R. McCarroll 


1. Introduction. In many problems of physics it is necessary to evaluate the 
spherical Bessel functions over a wide range of argument and up to a high order. 
For example, their evaluation is necessary in the solution of the integral equations 
of atomic scattering, as described in the work of Frazer [1]. 

One method is to generate the functions by means of recurrence formulas [2], [3], [4], 
which basically provide a large number of orders at a single value of the argument. 
A method which in essence provides a single order for a given range of argument, and 
which is suitable for use in automatic computations associated with atomic scatter- 
ing, is to expand the spherical Bessel functions in terms of Chebyshev polynomials. 


2. The Method. We introduce the shifted Chebyshev polynomial T,,*(z) which 
satisfies the following differential equation given by Lanczos [5]. 








_ 9 GT. (28-1) dTs® , amg 
(1) (z—z) 7A — hue +nT,* =0 0 


lA 
N 
lA 


The spherical Bessel function 


(2) jx) = (4/22)? Jrsin(2) 


is expanded in series of T,,*(z) for z S 1 as 


(3) jx) = (2/2)"N, D> AaT,.*(z) 
where 
(4) z= 2/p 


and where N, is the normalization factor given by 

(5) N, = x2 T(r + 3/2) >, (—1)"Aal, 
chosen so that as x — 0 

(6) jx) — w'?(x/2)*/2 V(r + 3/2). 


The parameter p is chosen according to the required z-range. 
Substituting the expression (3) into the differential equation satisfied by j-(z), 
namely, 


da 2d - t+) ; ast 
(7) E + = ae +1 mos Bias j(x) =0, 
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TaBLe 1 


Tables of Expansion Coefficients 


je(z) = (2/2)"N, Do, AnT'n*(z*/100) 


-10S323510 









































" a. | An 
measiad, tai 
| 
0 1.0000 0000 0000 1.0000 0000 0000 
1 —1.8370 2985 7106 —1.7335 5787 9572 
2 1.6181 7789 5224 1.7037 9637 8855 
3 —2.4956 2512 4954 —1.2117 2623 8820 
4 1.6832 2166 9128 0.4966 0881 5456 
5 —0.5891 1635 2032 —0.1261 1219 3614 
6 0.1275 7849 3735 0.0216 6172 6076 
7 —0.0189 6827 0697 —0.0026 8726 9369 
8 0.0020 6883 1574 0.0002 5255 1305 
9 —0.0001 7326 9614 —0.0000 1863 3794 
10 0.0000 1152 2733 0.0000 0110 9455 
11 | —0.0000 0062 4256 —0.0000 0005 4481 
12 | 0.0000 0002 8116 0.0000 0000 2246 
13 | —0.0000 0000 1070 —0.0000 0000 0079 
14 | 0.0000 0000 0035 0.0000 0000 0002 
15 —0.0000 0000 0001 
js 
n a An 
0 | 1.0000 0000 0000 1.0000 0000 0000 
ee tt er 
2 | . ; 
3 | —0.7040 0124 0486 —9.4322 1712 7394 
4 0.2158 5586 1952 0.1074 1449 1493 
5 —0.0438 6700 5066 —(}.0183 5466 7482 
6 0.0063 0718 4207 0.0022 8121 8406 
7 —0.0006 7512 0822 —).0002 1543 2710 
8 | 0.0000 5593 1115 0.0000 1599 1357 
9 —0.0000 0369 5921 —0.0000 0095 8213 
10 0.0000 0019 9498 0.0000 0004 7352 
ll —0.0000 0000 8968 —0. 9000 0000 1964 
12 0.0000 0000 0341 0.0000 0000 0069 
13 —0.0000 0000 0011 | —0.0000 0000 0002 
s i | An 
0 1.0000 0000 0000 | 1.0000 0000 0000 
1 —1.6836 0574 1497 —1.5488 1806 4538 
2 | 0.8722 0517 6005 0.6859 0242 4243 
3 —0.2782 8514 3648 —0.1872 6092 6662 
4 0.0585 4062 0216 0.0342 7629 1788 
5 —0.0086 6953 0055 —§, 777 
6 0.0009 5173 1760 0.0004 4233 6046 
7 —0.0000 8057 5832 —0.0000 3398 2101 
8 | 0.0000 0542 5189 0.0000 0209 5486 
9 —0.0600 0029 7662 —0.0000 0010 6104 
10 0.0000 0001 3573 | 0.0000 0000 4494 
11 | —0.0000 0000 0523 —0.0000 0000 0162 
12 0.0000 0000 0017 0.0000 0000 0005 
| ] 
n A. An 
0 | 1.0000 0000 0000 1.0000 0000 0000 
1 —1.4171 3242 2748 —1.2972 4767 3744 
2 0.5479 5705 4665 | 0.4451 8041 2800 
3 —0.1310 4234 0105 | —0.0948 5564 9030 
4 0.0212 7626 4218 0.0138 5756 5596 
5 —0.0025 0394 3447 —0.0014 8187 0216 
6 0.0002 2390 8476 0.0001 2144 5196 
7 —0.0000 1575 8616 —0.0000 0789 1023 
8 0.0000 0089 6941 0.0000 0041 7246 
9 —0.0000 0004 2186 —0.0000 6001 $329 
10 0.0000 0000 1669 0.0000 0000 0680 
11 | —0.0000 0000 0056 —0.0000 0000 0022 
12 


0.0000 0000 0002 








0.0000 0000 0001 
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Tables of Expansion Coefficients—Continued 
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0 1.0000 0000 0000 | 1.0000 0000 0000 
1 —1.1911 3494 9922 —1.0981 6702 5095 
2 0.3675 0203 3005 0.3078 0618 7533 
3 —0.0706 7640 7499 —0.0539 7859 5631 
4 0.0093 9388 4497 0.0065 8545 8180 
5 —0.0009 2119 2296 —0.0005 9665 6926 
6 0.0000 6972 0100 0.0000 4196 

7 —0.0000 0420 9178 —0.0000 0236 6972 
8 0.0000 0020 7889 0.0000 0010 9704 
9 —0.0000 0000 8569 —0.0000 0000 4260 
10 0.0000 0000 0300 0.0000 0000 0141 
ll —0.0000 0000 0009 —0.0000 0000 0004 
n - ym 

0 1.0000 0000 0000 1.0000 0000 0000 
1 —1.0168 7201 6138 —0.9456 3737 0019 
2 0.2611 5956 8765 0.2241 3285 9656 
3 —0.0421 0983 2180 —0.0334 5827 7279 
+ 0.0047 5003 6264 -0035 1079 6899 
5 —0.0004 0010 7956 —0.0002 7638 9359 
6 0.0000 2629 5870 -0000 1705 0581 
7 .0000 0085 0428 
8 ' 

9 
10 
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—0.0000 0139 1880 —0 
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0 
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0.0000 0006 0782 0000 0003 5115 
—0.0000 0000 2232 —0.0000 0000 1223 
0.0000 0000 0070 -0000 0000 0036 
—0.0000 0000 0002 —0.0000 0000 0001 
jiz jus 
An A 
1.0000 0000 1.0000 0000 0000 
—0.8829 7010 1098 —0.8275 7427 4817 
0.1943 1493 5452 0.1699 8535 
—0.0270 1121 0330 —0.0221 1284 8562 
0.0026 5020 1752 0.0020 3772 5111 
—0.0001 9588 4596 —0.0001 4196 6316 
0.0000 1138 8223 0.0000 0780 5415 
—0.0000 0053 7112 —0.0000 0034 9193 
0.0000 0002 1036 0.0000 0001 3008 
—0.0000 0000 0697 —0.0000 0000 0411 
0.0000 0000 0020 0.0000 0000 0011 
ju ju 
n An 
1.0000 0000 0000 | 1.0000 0000 0000 
—0.7783 5629 2813 —0.7344 0438 9543 
0.1498 9703 5766 | 0.1331 3186 9771 
—0.0183 2686 5868 —0.0153 5579 5806 
0.0015 9235 3331 0.0012 6227 9463 
—0.0001 0493 0348 —0.0000 7891 7018 
0.0000 0547 2793 0.0000 0391 5367 
—0.0000 0023 2885 —0.0000 0015 8871 
0.0000 0000 8272 0.0000 0000 5393 
—0.0000 0000 0250 —0.0000 0000 0156 
0.0000 0000 0006 


0.0006 0000 0004 
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TABLE 2 
Table of Normalization Factors 











r | Nr 

0 | 10-' -1.0670 1130 3957 
1 10- -1.0588 0224 7273 
2 10 -5.0977 3196 1602 
3 10 -1.7016 2862 5983 
4 10-* -4.3387 2970 2098 
5 10 -8.8939 6388 4585 
6 10-4 -1.5178 7598 0079 
7 10 -2.2135 3657 5215 
8 10* -2.8143 4241 8551 
9 10“? -3.1696 2427 4560 
10 10-* -3.2028 4879 3012 
11 10 -2.9343 5222 5384 
12 10-”-2.4587 5268 1163 
13 | 10-"- 1.8981 3134 3689 
14 | 10-”- 1.3584 9275 4731 
15 10-"-9.0623 7664 9006 





and using the properties of the shifted Chebyshev polynomials, we obtain a set of 
simultaneous linear algebraic equations which may be solved for the ratios A;/Ao, 
A2/Ao, --- . The solutions may be normalized, using the factors given by (5). 


3. Calculations and Results. The calculation of the coefficients A, has been 
carried out on a DEUCE digital computer for p = 100, and r = 0 to 15. For all 
values of r, Ao has been chosen as 1.0, and the ratios A,,/A» computed accordingly. 
Tables of the coefficients, together with the corresponding normalization factors are 
presented herein. The expansion coefficients are given to 12 decimal places to insure 
that for the range of x considered the spherical Bessel functions should be accurate 
to 10 significant figures. 

As can be seen from the tables, the convergence of the coefficients is very rapid; 
if the Chebyshev expansion and Taylor series are curtailed after n terms, the ratio 
of the (n + 1)th terms is about 1/2”. 

The authors wish to thank Mrs. N. Scott and Miss M. Wilson for assistance 
with the computations. 
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Gaussian Quadrature with Weight Function x" 
on the Interval (— 1, 1) 


By Harry A. Rothmann 


Most of the existent literature dealing with Gaussian quadrature formulas is 
based on the assumption that the relevant weight function is of constant sign in 
the interval of integration. This note deals with the weight function z” on the 
interval (—1, 1) and indicates the extent to which the existent theory can be 
generalized in that case. 


Gaussian quadrature formulas on the interval (—1, 1) have the form 
1 m 
(1) [ w(2)f(2) de = Do Wafer) + E 


where the weights and abscissas are chosen to ensure a degree of precision 2m — 1 
(i.e., exactness for all polynomials of degree not exceeding 2m — 1). 

One method of determining the abscissas z, involves obtaining a set of poly- 
nomials, ¢; , $2 , ¢3, -** , Such that each is orthogonal to all polynomials of inferior 
degree relative to the weight function w(x) over the interval (—1, 1). That is, for 
the mth-degree polynomial ¢,,(2), 


(2) [, w(2)on(2)an-a(2) de = 0 


where gm-i is an arbitrary polynomial of degree m — 1 or less. The m abscissas 
x, then are the zeros of the polynomial ¢,, . 
Upon defining 


(3) w(2)da(z) = £ al) 


the requirement that ¢,, be of degree m implies that U,,(z) must satisfy the dif- 
ferential equation 


(4) 0 | vets) ae 


dx™ | w(x) dx™ 
in the interval (—1, 1). From expression (2), the 2m boundary conditions 
(5) Um(+1) = Un'(+1) = +--+ = Un (+1) = 0 
can be obtained. When w(x) = x”, the solution of (4) is found to be of the form 
Un(xz) = x™w(x) [co + at +--+ + emt”) + do + da + --- + daar” 


from which ¢, can be obtained. It is convenient to impose the additional normalizing 
property 


(6) omn(1) = 1. 
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From the theory as given in [1], the weight corresponding to the abscissas :, 
for the N-point formula is 


abe o(e) & 
(7) Wise stay 1, z) : 


When the weight function w(z) is of constant sign, (7) reduces to 





Ay Yn-1 Aw41 YN 
8 W ‘4 = 7 = r 
(8) Nit” Ax-1 On (@x)on—1(2e) Aw ow (Xx) Ow4i(Ze) 


where Ay is the coefficient of x” in oy and 








(9) vy = Ay [ w(x)aoy(x) dz. 

The error function for the N-point formula can be expressed as 

(10) Ey = [ w(x)f [z1,21,°++ , tv, tw, 2)ey (x) de 

with my(x) = (x — %)(x% — a2) --- (@ — 2&y) where the z; are the zeros of gy . 
The expression (10) reduces in the case when w(x) is of constant sign to 

(11) Ey = sant 


where f(x) is assumed to have 2N continuous derivatives in (—1, 1) and 7 is some 
value in that interval. 


When w(x) = x”, the procedure outlined above can be employed to obtain 
the following 


Dd) = 1 ®,(zx) = 72 


$,(x) = 5 ((2n + 3)a” — (2n + 1)] 


(x) = 5 (Qn + 5)a* — (2n + 3)z] 


$,(7) = ap l2 n+ 7)(2n + 5)a* — 2(2n + 5) (Qn + 3)2” + (2n + 3)(2n + 1)] 


$;(z) ie 





—s 7)a° — 2(2n + 7)(2n + 5)2* + (2n+5)(2n + 3)a]. 


From these, the forms 


N N 
(12) tan) = sete () (—1)' [1 an +200 ~ 1) + 25 - at 


t 


(m = 2N) 


(13) tale) = ts (")(—v'T ba + 20W 1) +9) 412" 


=~ t=O 





(m = 2N +1) 
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can be deduced [2] by an inductive procedure and shown to satisfy the desired 
properties (2) and (6). Since the weight function is of constant sign in the interval 
(—1, 1), a known theorem [1, page 171-2] states that the zeros of the polynomials 
are all real, distinct, and lie in the interval (—1, 1). 

To display the formulas for the weights and error, the values Ay and yy must 
be obtained. From (12) and (13), it is easily seen that 


Aw = mie Hosen 1) 
(14) 
Awwsi = _ a | I (2n + 2N + 2j +1). 
The value 
2 
(15) ip ee 


2n + 2N +- 1 
can be established by an inductive proof [2]. Thus, the weights are now expressed 
in the form 


1 —2 
Wawa = N@ oy (2) P2y_1( 2) m (2n +2N + 1) @oy( a) Pon 4s (24) 
Wows ts* a ons , tas, — > 
wor (Qn + 2N + 1)@on41 (2p) Poy (ay) CN + 1) @oy41( rp) Boys2(x) 








Since w(x) is of constant sign, the error is given by (11) which can be written 
in the form 


age _2(N12" YF i Ae 
(4N)!(2n + 4N + 1) NH (2n + 2N + 2j — 1)? 


2(N12")' FO"? (n) 





Fonsi = " - ~~ ———_ 
(4n +2)! (Qn 4+ 4N Pe FE (-1<<1). 
j=l 


2n+l1 


When w(x) = 2x””, it is found [2] that the formulas for an odd number of 
abscissas do not exist, but that the polynomials of even degree satisfying properties 
(2) and (6) do exist, and have the form 


(16) 62x (x) = weds \-'D I] [2n + 2(N - 1) + 2 + 1) 2°" . 


j=l 


Since x6ey(x) = en4;(x), the zeros of Oey are all real and distinct and lie in the 
interval (—1, 1), and, in fact, are the zeros of 42y,; when the zero x = 0 is sup- 
pressed. 


The weights for the 2N-point forntula involve the polynomial 4.,,; in the 


following manner 
1 
ont+1 Goy( 2) 
Weve = L [ pe ds 


ow (ay) I—- X% 





166 HARRY ROTHMANN 


1 / 
1 on Pon 4i1\ } 
SS & 
1 








~ Oaw( ay) t— 2 


which reduces to 


- — Aanse Yen+1 = Aoy+1 Yan 
(17) Won.x _ 7 am 7 
Aow+iOon(t,)Pewye(Ze)  AanOow( x.) Pow( re) 


where the A’s and 7’s refer to the polynomial ® and are given by (14) and (15). 
After the substitution of (14) and (15), (17) reduces to 





ot 2 


(N + 1) Ooy(2,) Pew 42( Xe) a (Qn + 2N +} 1) 6224) ey( 24) ; 
From (10), the error term E2y can be written as 


Won.x = 
































1 
E2y = [ ati [ay ,1,°** , Xan, Len , ZI] wn (2) dx 
TABLE 1 
Gaussian Quadrature with Weight x" 
m n Weights Abscissas Error Coefficients 
| |— —— 
2 0 | — 1.0000000 | + .5773503 7.4 x 10° 
. 4 . 3333333 | + .7745967 1.9 x 10° 
2 2000000 | + 18451543 7.6 x 107 
3 "1428571 | 4 .8819171 3.7 x 10+ 
4 1111111 | + .9045340 2.1 x 10+ 
5 | 09090909 | + .9198662 1.3 x 107 
3 0 | 8888889 0 6.3 x 10° 
5555556 | + .7745967 
1 . 1066667 0 | 2.5 x 10° 
. 2800000 + .8451543 
2 .03265306 0 1.2 x 10° 
. 1836735 + .8819171 
. 01410935 0 7.1 x 10-¢ 
. 1358025 + .9045340 
4 007346189 0 4.4 x 10-6 
. 1074380 + .9198662 
5 004303389 0 2.9 x 10-¢ 
08875740 + .9309493 
4 0 . 3478548 + .8611363 2.9 x 107 
6521452 + .3399810 | 
Rated 1945553 + .9061798 7.3 x 10° 
. 1387780 + .5384693 
2 | . 1343622 + .9290483 | 2.5 x 10-8 
| 06563784 + .6399973 
4 . 1024498 | + .9429254 1.0 x 10° 
| .04040730 | + .7039226 
= a .08273203 | + .9522526 4.9 x 10° 
| | .02837808 | + .7482524 
| 5 | (06935661 | + .9589554 2.6 x 107 
+ 





-02155248 


- 7809074 

















GA 


wl 
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which reduces after one integration by parts to 

1 z 

Ey = —[ sles 771 pr MT > Ten » Zan , 2, 2) [ ogy (2) dz dz. 
1 
Now consider the function 
A(z) = [ a giy (x) dr. 
1 


Since A’(x) is negative for z < 0 and positive for z > 0 and A(+1) = 0, it follows 
that A(z) is of constant sign in (—1, 1) and hence 





1 z 
E,y = —flti,--- , tev, €, 8 [ [ sine) dx (-1<é< 1) 
(18) ‘ 
fn) 1 a 
= (GN + 1)! 2 Kow4i(x) dx (-1 <9 <1) 


where f(z) is assumed to have 4N + 1 continuous derivatives in (—1, 1). 
The function Key; is the monic polynomial consisting of the linear factors of 
.y4; and therefore (18) reduces to 





Eu: = for (9) ren 
(4N + 1)! Abwas 


i.e. 2(N12")*f“"*" (9) 
Eon = 





N ° 
(4N + 1)1(2n + 4N + 3) [J (2n + 2N + 27 + 1)’ 


j=1 


TABLE 2 
Gaussian Quadrature with Weight x2" 














m n Weights Abscissas | Error coefficients 
2 | 0 + .4303315 + .7745967 3.8 x 107 
1 + .2366432 + .8451543 | 1.5 x 10° 
| 2 + .1619848 + .8819171 | 7.5 x 10° 
3 + . 1228380 + .9045340 | 4.2 x 10° 
4 + .09882860 | + .9198662 2.6 x 10° 
| 5 + .08262864 | + .9309493 1.7 x 10° 
| 
4 | 0 + .2577268 + .5384693 8.1 x 10° 
+ .2146984 + .9061798 

1 + . 1025596 + .6399973 2.8 x 10° 

| + .1446234 + .9290483 
2 + .05740305 + .7039226 1.1 x 10° 

+ .1086511 + .9429254 
| 3 + .03792714 + .7482524 5.5 x 107° 

+ .08688035 + .9522526 
4 + .02759928 + .7809074 2.9 x 10> 

+ .07232517 + .9589554 
5 + .02137648 + . 8060023 1.6 x 107° 

+ + .9640060 


.06192242 
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In Tables 1 and 2, the weights and abscissas are given for m = 2, 3, 4, for w(x) = 
x” and for m = 2, 4, in the case w(x) = 2°"*". The values for the weights and 
abscissas are correct to the 7 figures given. The coefficient of f() in the error 
term is also given. These error coefficients are correct to the two figures given. 
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A Table of Generalized Circular Error 
By Harry Weingarten and A. R. Di Donato 


1. Introduction. This note provides an abbreviated table (Table 1) giving sclu- 
tions for the value of K satisfying 


1 ifx ,¥ . 
sama, I, few[ - $(B + 4) aca che 


where R is the circle z* + 7 = K’s,’, 0, = o, and c = o,/c,. K has been com- 
puted for c = 0(.01)1 and P = 0(.01).99. The table provided here will not contain 
all the results because of space limitations. The complete table is available upon 
request directed to either author. It differs from the extensive one in [1] which also 
gives numerous applications and a wide bibliography of the bivariate normal 
distribution. 





2. Application. When P = .5 and c = 1 we obtain the cPE (circular probable 
error) relationship used in ballistic studies. In this case K = 1.17741, which may 
easily be found without the table in this note. When c ¥ 1, however, (which is the 
usual case) it is still of interest to find the circles within which impacts will occur 
with given probabilities, rather than the ellipses. For any particular P and c the 
value of K in the table multiplied by o, is the required radius. 


3. Statistical Interpretation. This kind of problem has been widely considered 
as indicated by the references in [2], where the approach is differently oriented, 
being concerned with the general problem of the distribution of quadratic forms. 
Essentially we consider here the cumulative distribution in tabular form of the 
random variable, 

Z=X*+Y 


where X and Y are independently and normally distributed with zero means and 
variances o, and a, . (Z does not, of course, have a x’ distribution unless ¢, = oy = 
1.) In [3], Chapter 27, there will be found application of such results to the specifica- 
tion of regions of type C in the testing of hypotheses. 


4. Analysis. This section will detail the computational and numerical analysis 
aspects of the preparation of the table. 
The probability integral under consideration is given by the following equation: 


1 1f/x , 
(1) P(K, ¢.,¢y) = — [Joo] -3(% +¥) lac dy 
aad z y = y 


where the region, R, is specified as a circle with its center at the origin and with 
radius Ko, . The use of polar coordinates transforms (1) to 


Qn x 2 2 
(2) P(K,c) = = [ l exp | — 40 (Ads ae = cos 20) |» dp dé 
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TABLE 1 
The Generalized Circular Probable Error K 


AND A. R. DI DONATO 








05 


-10 


15 





















































K 25 30 35 40 AS 50 
-05 | 0.08149) 0.10697) 0.12806 0.16251) 0.17730) 0.19097) 0.20375) 0.21579) 0.22721 
-10 | 0.13631) 0.16328) 0.19017 0.23662) 0.25701; 0.27599) 0.29383| 0.31070) 0.32675 
-15 | 0.19590) 0.21757) 0.24565 0.29897) 0.32313) 0.34585) 0.36734! 0.38777 0.40729 
-20 | 0.25834| 0.27454) 0.30048 0.35690) 0.38367) 0.40917) 0.43349] 0.45676) 0.47910 
-25 | 0.32259] 0.33506) 0.35715 0.41348) 0.44188) 0.46941) 0.49596] 0.52165) 0.54624 
.30 | 0.38858) 0.39867) 0.41685 0.47050) 0.49965) 0.52853) 0.55677) 0.58424) 0.61093 
.35 | 0.45653) 0.46500) 0.48004 0.52924| 0.55829) 0.58788) 0.61731| 0.64626) 0.67463 
-40 | 0.52679) 0.53409) 0.54679 0.59073) 0.61889) 0.64854) 0.67866) 0.70872) 0.73846 
-45 | 0.59986) 0.60623) 0.61721 0.65585) 0.68244) 0.71154) 0.74184! 0.77260) 0.80339 
.50 | 0.67635) 0.68199) 0.69163 0.72543) 0.74994| 0.77788) 0.80785) 0.83890) 0.87042 
.55 | 0.75707| 0.76210) 0.77066 0.80039} 0.82243) 0.84870) 0.87782| 0.90870) 0.94060 
-60 | 0.84311) 0.84761) 0.85527 0.88142) 0.90113) 0.92532) 0.95307) 0.98332) 1.01520 
.65 | 0.93593) 0.93998) 0.94685 0.97008) 0.98751) 1.00939) 1.03532] 1.06444) 1.09586 
-70 | 1.03764) 1.04129} 1.04748) 1. 1.06822} 1.08361) 1.10311) 1.12685) 1.15433) 1.18481 
.75 | 1.15144] 1.15473) 1.16029) 1. 1.17884} 1.19246} 1.20968) 1.23100) 1.25637) 1.28534 
.80 | 1.28253) 1.28548) 1.29046) 1. 1.30704) 1.31908) 1.33421) 1.35302) 1.37588] 1.40275 
.85 | 1.44040) 1.44303) 1.44746) 1. 1.46215) 1.47277| 1.48599) 1.50233) 1.52238) 1.54653 
.90 | 1.64561) 1.64791) 1.65179) 1. 1.66461} 1.67383) 1.68523) 1.69918) 1.71626) 1.73708 
.95 | 1.96060) 1.96253) 1.96578) 1. 1.97651| 1.98420) 1.99366) 2.00514) 2.01902) 2.03586 
.96 | 2.05436) 2.05620) 2.05930) 2. 2.06953) 2.07686) 2.08587) 2.09679] 2.10995) 2.12588 
.97 | 2.17067| 2.17241) 2.17534) 2. 2.18502} 2.19195) 2.20045) 2.21075] 2.22314) 2.23806 
.98 | 2.32689) 2.32851) 2.33124 2. 2.34026) 2.34672) 2.35464) 2.36421) 2.37569| 2.38948 
.99 | 2.57632| 2.57778) 2.58025) 2. 2.58839) 2.59421 ers 2.60995) 2.62025) 2.63257 

KK 55 60 65 .75 80 85 90 95 1.0 
.05 | 0.23810) 0.24852) 0. 0. 0.27753) 0.28657| 0.29534) 0.30388) 0.31219) 0.32029 
.10 | 0.34210) 0.35683) 0.37101) 0. 0.39798) 0.41085) 0.42336) 0.43555] 0.44744) 0.45904 
.15 | 0.42601) 0.44402) 0.46142) 0. 0.49458) 0.51045) 0.52591) 0.54099) 0.55571) 0.57012 
-20 | 0.50060) 0.52136) 0.54145) 0. 0.57990) 0.59835) 0.61636) 0.63396] 0.65118) 0.66805 
.25 | 0.57012) 0.59326) 0.61573) 0. 0.65888) 0.67967) 0.69999) 0.71989] 0.73939) 0.75853 
.30 | 0.63688) 0.66213) 0.68672) 0. 0.73418) 0.75712) 0.77960) 0.80166) 0.82331) 0.84460 
.35 | 0.70237) 0.72950) 0.75604) 0. 0.80748) 0.83246) 0.85699) 0.88110) 0.90483) 0.92821 
-40 | 0.76775) 0.79655) 0.82486) 0. 0.88004) 0.90696) 0.93346) 0.95958) 0.98534) 1.01077 
.45 | 0.83399) 0.86428) 0.89421) 0. 0.95291) 0.98170) 1.01013 1.03822) 1.06599) 1.09347 
.50 | 0.90207) 0.93365) 0.96505) 0. 1.02709} 1.05769} 1.08801) 1.11807) 1.14786) 1.17741 
.55 | 0.97303) 1.00569) 1.03841) 1. 1.10361} 1.13599} 1.16819) 1.20021) 1.23206) 1.26373 
-60 | 1.04810} 1.08162) 1.11549) 1. 1.18366} 1.21779) 1.25187) 1.28590) 1.31985) 1.35373 
.65 | 1.12888) 1.16298} 1.19781) 1. 1.26875} 1.30460) 1.34058) 1.37666) 1.41281| 1.44901 
.70 | 1.21752) 1.25187| 1.28742) 1. 1.36090) 1.39845) 1.43637| 1.47459) 1.51306) 1.55176 
.75 | 1.31724) 1.35143] 1.38739) 1. 1.46309) 1.50231) 1.54222) 1.58271) 1.62369) 1.66511 
-80 | 1.43320) 1.46668) 1.50262) 1. 1.58010} 1.62096) 1.66294) 1.70586) 1.74962) 1.79412 
.85 | 1.57477) 1.60677) 1.64206) 1. 1.72059} 1.76302) 1.80717) 1.85280) 1.89974! 1.94788 
.90 | 1.76212) 1.79152) 1.82511) 1. 1.90335) 1.94716) 1.99359) 2.04236 2.09321) 2.14597 
-95 | 2.05638) 2.08130) 2.11111) 2. 2.18580} 2.23029) 2.27908) 2.33180] 2.38812) 2.44775 
.96 | 2.14527) 2.16891! 2.19748 2.27054| 2.31491) 2.36413] 2.41782) 2.47565) 2.53727 
.97 | 2.25619) 2.27835) 2.30537 2.37617| 2.42021) 2.46978) 2.52455] 2.58415) 2.64823 
.98 | 2.40614) 2.42650) 2.45153 2.51895) 2.56226) 2.61202) 2.66799] 2.72983! 2.79715 
-99 | 2.64736) 2.66533) 2.68750 2.74916) 2.79069) 2.84010) 2.89743 


2.96249) 3.03485 
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where p, @ are the usual polar coordinates stretched by a factor o, and where 


(3) 0<%=c¢ a 


oz 


lA 
A 


Simple transformations reduce P(K, c) to 





1 2/2 7 
(4) P(K,c) = Lf mad | eA? 19 dw 
where 
i-¢ 1 ? 
(5) A =->* ond B= *tS. 


The integral over 6 in (4) is referred to in [4, (p. 46)], and may be expressed as 


(6) [ eve? dg = x1y( Aw) 


where J)(z) is defined by the following Taylor and asymptotic expansions re- 
spectively 


e) 1 2 2n 
(7) I(x) = (2) (5) 

e XS [(2n)P _, 
(8) I,(x) siinaal </2Qxx » 25"(n!)3 vt: 


The relations (7) and (8) are given in [4] on pages 20 and 58, respectively. 
Two computation schemes were used for computing P. if AK’ < 40 (an arbi- 
trary choice), then the following recurrence relation was used to compute P: 


2n —1(AY 1 ae) _pxtp (AK? (4) 2) 
(9) Tx = —- (3) T n-2 zl ( rn e = +n B)i Nai 








where 
1 —BK? /2 
(10) Ty = Be (l-—e ) 
and 
N 
(11) P= p> 7 AK’ < 40. 


If AK’ > 40 then the following recurrence relation was used to compute P: 


1 1 2 [(2n) !}" A 2\—@n-Dl2 (—K2/2 
(12) Mass = 375° Vein — 12 AK) * 
L»—l 
A 2n 


M 2n—1 


« 


~ 


t 


where 


(13) M, == Ra a “ =f em dw 
7 KiV 
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and 
me 

(14) P=1— DY May. 
n=O 


Equation (9) is obtained by substituting (7) into (4), transforming the upper 
limit on the W integration from K*/2 to AK’*/4, interchanging summation and 
integration, and then performing two successive integrations by parts on the integral 
that occurs as part of the general nth term of the series. The upper limit, N, of the 
sum that appears in (11) is determined when 


W | 
(15) | Taw| S€| D0 Ton | 
n=0 } 

where ¢ is chosen to the order of accuracy to which P is desired. 

The recurrence relation given by (12) is derived by substituting (8) into (4), 
interchanging summation and integration, and by considering the integral from 
AK’ to infinity rather than from 0 to AK’. Two integrations by parts on the integral 
that occurs as part of the general nth term of the series yield (12). The integer N’ 
is determined such that 


N’ | 
(16) | Mowr4i| Se » Mon+1 








The restriction of (12) to the region AK* > 40 insures at least eight-digit accuracy 
in P before the terms M.2,.; eventually begin to increase in magnitude. The e’s in 
(15) and (16) wer. vet at 10°. 

Inasmuch as equal intervals in P and c were desired, a Newton-Raphson pro- 
cedure was used to determine K for a given P and c; accordingly 


1 pons! 
“ [ e °"I,(Aw) dw — P 
0 


(17) Ky sas | +i S 





~ 
a —(BK? /2 AK, 
=k. ot Me 5 
Cc 





- 


where K,, represents the nth iterate for K. 

The efficiency and accuracy of the computation are indicated by the fact that 
the average time required to evaluate a K to eight significant digits for a given P 
and c was 50 milliseconds on NORC. The accuracy of the results was checked by 
evaluating the same K by both (9) and (12) in the region of AK” = 40. Thirty 
terms were used for this purpose. This region is where both series (11) and (14) 
require the largest number of terms, and consequently where truncation and round- 
ing errors should be the largest. Some further checks to insure eight-digit accuracy 
were obtained by evaluating some of the integrals by the direct application of 
Simpson’s Rule. The entire table presented herein required less than 30 seconds of 
computing time on NORC. 
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Polynomial Approximations to 
Integral Transforms 


“> 


By Jet Wimp 


1. Introduction. The symmetric Jacobi polynomials P{*'*’(z), orthogonal on 
the interval —1 < z < 1, are widely used for approximating functions, but the 
integral which defines the coefficients for the expansion of a function g(z) in these 
polynomials usually is quite difficult to evaluate. The problem is simplified if g(z) 
is an integral transform of the Fourier or Laplace type, since the kernel of the trans- 
form generates a series of the above polynomials. The coefficients in such cases are 
found to be Hankel transforms, which are widely tabulated. 

Examples include Chebyshev polynomial expansions of 1/(z + a)*, ¥(z + a), 
log T(x + a), Ci(x) and Si(z). 


2. Formulas When g(x) is a Laplace or Fourier Transform. The symmetric 
Jacobi polynomials [1, v. 2, p. 168] may be defined by 
(1) P(x) = ("3") Fil—n, n + 2a + 1; a + 1; 3 — 4a). 


A function g(x) satisfying certain conditions has the expansion 








(2) g(z) = 2d, AaPa* “(z), -l1s7rK1, 
where 

/ 1 
a) 4.2 SS et f olen — PPL) ae. 


PHT (n +a + IP ‘? 
Suppose now that g(x) is the Laplace transform of some f(t), 
(4) gz) = ef) = [ e“Pat = DY A, PS (2). 
0 n= 
To determine the A,’s replace the kernel of the Laplace transform by its Neumann 


series [1, v. 2: p. 98, No. (1); p. 175, No. (16); p. 174, No. (6); and the duplication 
formula for the gamma function]. 





—zr . n Ey a+) / (t) (a,a 
(5) é : se » (—) 2, aa Ps (x), 
(6) go a 2 (n + a + 4)T(n + 2a + 1) 
" Tin + a+ 1) : 


Then (4) yields 


_ _@e-1) [esi f(t) 
(7) A, =e nS 1 ” , 
\ y=on+a+1/2 
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(8) se(F()} = [- POIA yt) yt)" at 
%x{F(t)} denotes the Hankel transform of F(t) [2]. 
When a = —4}, (7) furnishes the coefficients for the Chebyshev expansion 
(9) g(x) = [ e*f(t) dt = Eo. T,(z), -l1s21, 
where 


If we replace ¢ by it in (5), we find that the same method is applicable when 
g(x) is a Fourier transform of f(t). We omit details, but the key results for the sine 
and cosine transforms are as follows. 


g(x) = a yon = y S, Ae.a)/ _— 

(11) ot) = b I cog @) tt = Do Pe (a), 1szs1 
where 

0, n even, 

2 S, = n—1) [ri t 
(12) ern mange : n odd, 
v=n+a+1/2 
and 
0, n odd, 
(13) C. = omi0, 0 FON ; n even. 
v=on+a+1/2 


3. The Chebyshev Expansion for 1/(y + a)". Let g(x) = [: 
Then 


to] + 


—k 
= +a] ‘ 


k 


(14) £"t9(2)} = Ep 


geome én f(t). 








Use (10) and let y = : * I Then T,(2y — 1) = T,*(y),0 S y & 1, is the shifted 
Chebyshev polynomial [3] and 
1 fS a(—)"(k+n— 1)! 5-5 
quits an P;- 
wtar~ \% (k= 1)! = 
~~ 2a + 1 ) 
a é 2 k/2 
- | —_——— | T,* / 0sy 81, > 0, 
|. r -| (wf (a + a) y a 


where P,”(x) is the Legendre function [1, v. 1. p. 120]. For k = 1, (15) agrees 
with a result of Luke [4]. 





SS a 






































Com we et OS Be! SS 
ne tee E ; 2 ers ee ek ct kl : | ' 
—- | = — —— Dies | | 10000000° — 
~ - = | _ _ | “ 10000000°— | #0000000" 
-_ + on | _ | on | 20000000" | 91000000 * — 
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4. The Psi and Log Gamma Functions. These examples show how a property 
of the Laplace transform may be used to advantage when applying (4) and (8). We 
know that 


(16) Le'f(t)} = g(x + a). 


If g(2) cannot be expanded in symmetric Jacobi polynomials, a in (16) can often 
be chosen so that g(z + a) has a convergent expansion. Let 


(17) g(x) = ¥™ (xz) = D™* log T(z). 


Since y‘”’ (x) has poles at zero and the negative integers, we cannot expand the 
function over —1 S z = 1. However, if 


(18) g(x) = y'"(x + a), 
then 
(19) f(t) = £"{g(z)} = (-)" eh — eT", 


and if Re(a) > 1, (7), and in particular (10), may be used since (18) is analytic 
for | z| < 1. Substituting (19) in (10) and expanding (1 — e~‘)~™ by the binomial 
theorem, we have 





























(20) C. = -. = k vf o1- 9 | 
k=o dx™ V2 —-—1 lpmkto 
Setting m equal to zero, we get 
af - 
(21) ~ [ (k +a)? — 1 « + a)” n21 
V(k + a? - 
TABLE 2 
Coefficients for the Series 
Ci(x) = / cs a = log(z) + >> AonTen (2), 0 < z<a 
r) n=0 
Si(z) = [ = tat = po» BonsrT on+1 (2) ; —@ = Zz = a 
0 t n=0 a 
: i 5 ra a=2 a=5 
: -decsaieaenineel vm stall peayiacensanneeipeeennnacune in 
Ain | Bony Ain | Bans 
0 0.13529 ) 62627 1.69809 09708 |—0.96313 15550 | 2.08578 21107 
1 | —.42327 51922 | —.09558 49521 |—1.13103 16550 | — .67042 59749 
S-.4 .01822 27219 .00295 78196 | -34661 70891 | .15186 68742 
3 — .00041 57650 | — .00005 14215 | — .05698 43620 | — .01861 43512 
4 | .00000 56716 | .00000 05642 .00537 47844 | .00138 96747 
5 — .00000 00511 | — .00000 00042 | — .00032 52237 | — .00006 95137 
6 |  .00000 00003 | a .00001 36729 |  .00000 24908 
7 — - | — .00000 04226 | — .00000 00671 
8 


| sink | os -00000 00100 |  .00000 00014 
| os | - — 00000 00002 on 
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If n = 0, (21) diverges, and for n = 1 the series is slowly convergent, but since 
T,(1) = 1, T,(—1) = (—)", we may solve for Cy and C;, in terms of higher com- 
putable coefficients, i.e., 


je, = MAE 4 Me 9): Fo, | 








(22) : E 
lo Wet —We-)_ Fo, 


Integration of the series defined by (21) yields a Chebyshev expansion for 
In T(x + a) because [3] 


nti(x) T,-(z) 
(23) [ tae) ae = $[ Pon) _ oan | +o. 


In Table 1 are listed coefficients for the Chebyshev expansions of ¥(x + a) and 
log T(x + a), a = 2(1)5, n = 0(1)15 to 8D. 





5. The Sine and Cosine Integrals. For examples of (11)—(13) let 
(24) g(x) _ (1 — cos ar)/z -[ f(t) sin xt 


g(z)  sinaz/x cos at? 
A. G<2< 6, 

(25) f(t) = 
0, sxis < @. 


Using [2, v. 2, p. 333, No. (1)] to evaluate (12) and (13) for a = —3, we find that 


0, n even, 
(26) Be =) piesa F 
[te 1) (ri/2] & J n+on41(@), n odd, 


0, n odd, 
(27) C. ss cae ae 
Slr FS Svnsle), n even. | 
k=0 


Let a = 2 and 5 in (26) and (27), and use [1, v. 2, p. 145, No. (6)] and (23) to 
obtain the expansion whose coefficients are listed in Table 2. 

The author wishes to thank Yudell L. Luke and also Charles Himmelberg and 
Jerry Fields for their helpful suggestions during the preparation of this paper. 
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On Stability Criteria of Explicit Difference 
Schemes for Certain Heat Conduction Problems 
with Uncommon Boundary Conditions 


By Arnold N. Lowan 


Abstract. Stability criteria are derived for the explicit difference schemes appro- 
priate to the following problems: 1) heat conduction in a slab in contact with a 
well stirred liquid; 2) heat conduction in a slab radiating to one face of a thin slab 
with infinite thermal conductivity, the other face of which radiates into a medium 
at prescribed temperature; 3) heat conduction in a cylinder radiating to the inner 
surface of a thin coaxial cylindrical shell with infinite thermal conductivity, the 
outer surface of which radiates into a medium at prescribed temperature. 


Although the exact analytical solutions of certain problems in heat conduction in- 
volving complicated boundary conditions are known, the complexity of the analyt- 
ical expressions is often such as to make them impractical for the numerical evalu- 
ation of the solutions. This is for instance the case of the writer’s solution of the 
problem of “heat conduction in a solid in contact with a well stirred liquid” [1]. It 
is also the case of the problems treated by Walter P. Reid and dealing with the heat 
conduction in a semi-infinite solid (or cylinder) when the boundary surface radiates 
to one boundary surface of a thin slab (or thin cylindrical shell), with infinite thermal 
conductivity, the other boundary surface of which radiates into a medium at pre- 
scribed temperature [2], [3]. 

To obtain numerical answers to the above problems it is expedient to evaluate 
the solutions of the appropriate explicit difference analogs. The object of this re- 
port is to derive the stability criteria for the difference schemes appropriate to the 
problems above-mentioned. 

Consider first the problem of heat conduction in a slab one face of which is in 
contact with a well-stirred liquid. For the sake of concreteness we shall first assume 
that the other face is kept at 0°C. The mathematical formulation of the problem is 
as follows: 


aT aT 
SS. eo <r< 
(1) at * ag er ga ptt 
(2) T(x, 0) = f(z) 
(3A) T(0,t) =0 
oT c,d\ aT oT 
(4) = = (2e¢) =° "a (say) for z =a. 


In (4) we have the boundary conditign appropriate to the case where the face 
x = ais in contact with a layer of a well-stirred liquid of width d, density pp and 
specific heat co. The constants K and k are the thermal conductivity and thermal 
diffusivity of the slab, respectively. 
Received August 17, 1959; revised November 4, 1960. 
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We shall consider the problem of stability of the explicit difference scheme 
(5) Taoan = TTaan t+ (1 — 2r)Tan trTasin; m=1,2,3,---M 


where Tn. = T(mAt, nAt);r = kAt/(Azx)’, and a = (M + 1) Az. The difference 
analog of (4) is 





‘ T uun ws T w+i.n al T w+i.n+1 et T uttin 
(6) ee oe At 
whence 
At 
(7) Tusine = 8Tun + (1 — 8)T tin; = ——. 
oAr 


Since At = r(Ax)*/k = r(Ax)*pc/K where p and c are the density and specific 
heat of the slab, it follows that s = r(Ax/d)(pc/poco). Since usually Ar < d and 
c < ¢ it is reasonable to assume that Az/d-pe/poco < 1 and a fortiori s < 1, if we 
anticipate the fact that the criterion for the stability of the explicit difference 
scheme under consideration is r < 3. 

The system of M + 1 equations consisting of (5) and (7) may be written in 
the matrix-vector form 


(8) Tau = A 7; 


where T, and T,,,; are the (M + 1) dimensional vectors whose components are 
the temperatures T,,,,, and T'm.41 , respectively, with m = 1, 2,3, --- M+ 1 and 
where 


(9) gg ye Bey te LOE pe ee 








8 1—s| 


Examination of the matrix A shows that ifr < 4 so that 1 — 2r = 0, then the 
largest of the sums of the absolute values of the elements of each of the first M 
rows of A is equal to 1. Furthermore, we have seen that it; is reasonable to assume 
that s S 1 or 1 — s 2 0; accordingly, the sum of the absolute values of the ele- 
ments of the last row of A is again = 1. The condition for the stability of the dif- 
ference scheme under consideration is thus satisfied [4]. 

In the above treatment we assumed that the face x = ( is kept at 0°C. If instead 
the temperature at x = 0 is prescribed, so that 


(3B) T(0, t) = $(t) 


it is readily seen that the stability criterion r < 4 we found above is valid, since 
the error vector E,, satisfies the difference equation (8) and since it vanishes on 
= 0. We define the error vector E, as the vector T, — T,*, where T, is the 


“true” solution of (8) and T,,* is the vector generated by successive application of 











—_— OS Oo 


n 
e 
yf 
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(8) when during some time step an inaccurate vector T,* (where k can of course 
be = 0) is used in lieu of the true T, . Ii may be briefly mentioned that the above 
stability criterion is equally valid in case of boundary conditions of the form 


oT 


(3C) —=0 for z=0 
Ox 

(3D) _.. hT for xz =0. 
Ox 


In the case of the boundary condition (3C) the first element of the first row of the 
matrix A is replaced by 1 — r; in the case of the boundary condition (3D) the 
element in question is replaced by 1 — 2r + ar where a = 1/(1 + hAz). In either 
case the largest of the sums of the absolute values of the elements of the rows of 
A is still equal to one and therefore the stability condition is satisfied ifr < }. 
For a discussion of convergence the reader is referred to [4], Section VI. 

Consider now the problem of heat conduction in a slab subject to the conditions 
stipulated in the first of Reid’s articles mentioned above. The mathematical formu- 
lation of the problem is as follows: 





aT eT 
(10) % = *s 0szs4a,t>0 
(11) T(x, 0) = f(x) 
(12) T(0,t) =0 
(13) -K(27) =m [T(a,t) — v(t) 

Ox z=a 

(14) po cod <? = hy [T(a,t) — v(t)] — hav(t) 
(15) v(0) = V. 
In the above: 


K = thermal conductivity of the thick slab 
k = thermal diffusivity of the thick slab 
h, = coefficient of heat transfer between the two slabs 
he = coefficient of heat transfer between the thin slab and the surrounding 
medium (with temperature zero) 
v(t) = temperature of thin slab 


po , Co and d are the density, specific heat and width of the thin slab. 
If we put 





. hy _ ‘ hp _ : Poco d _ 
(16) K =» .s K =d 
equations (13) and (14) assume the form _ 

(13*) sa (37) = b[T(a, t) — v(t)] 
Ox Jama 
(14*) oe = UT (a,t) — v(t)] — cot). 


ot 
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We shall investigate the stability of the explicit difference scheme 


Toot = Tae Fp... — 2T ae + Tass) 





At ~ (Ax)? 
or 
(17) eee = To. + (1 _ 2r)T mn -$- TT att.2 ’ = = 1, 2, 3, --- M 
where 
kAt a 
T m,n = T(mAz, ndt), r= (an)? and Ar = Mq1° 
The difference analogs of (13*) and (14*) are 
T i ag T n 
py one = BCT 41.0 ial Un) 
or 
1 bAr . 
(18) T u+in = 1+ baz Tun + i+ bAr Vn = PT un + Qn (say) 
and 
a “ = bT usin — (b +0) 
or 


_ (b+) At 


Co 


bAt 
(19) ty = — Putin + [1 |e. = aT yiin + Brn (say) 


where we have written v, and v,4; for v(nAt) and v[(n + 1) Ad], respectively. If from 
(18) and (19) we eliminate T y4;,. we get 


(20) Va41 = Aptun + (B + ag)un. 
If we rewrite (18) in the form 
(18*) T w+i.n+i = PT venti + QMn+i 


and eliminate v, and vp; from (18*), (18), and (19), we get 
(21) T w4inti = (8 + aq) T 41,0 + PT uin4a a PBT un . 


If in (21) we replace T' y,.4; by its expression from (17) with h = M, we ultimately 
get 


Pusan = PT wan + p(l — 27 — B)T un + (pr + B+ ag)T urin 
= Fa tite + QT un + ST wisn (say). 


(22) 


Starting with the values of T’,,,, for nm = 0 equations (17) and (22) will generate in 
succession the values of 7'»,, forn = 1, 2, --- and m = 1, 2,3, --- M + 1. Simi- 
larly, starting with v(t) = V for ¢ = 0 equation (19) will generate in succession the 
values of v, = v(ndAt). 
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The system of M + 1 equations consisting of (17) and (22) with the errors E,.., 
written for the T,,,,’s may be written in the matrix-vector form 
(23) E41 = AE, 


where E,, and E,,,; are the (7 + 1) dimensional vectors whose components are the 
errors in the values of the temperature at the M + 1 lattice points mAz, m = 
1, 2,3, --- M + 1 at the times nAt and (n + 1) At and where 








Ti—2r r ‘ 
r 1-—2r fr 
(24) ae r 1-—2r fr 
r i-2? ff 
' Pe si 
If we assume 
(25) 0<rs} 


it is clear from (24) that the sum of the absolute values of the elements of the first 
M rows of A is equal to 1. If in addition we assume that 


(26) 1-—2r-—620, pr+B6B+aq20 
then, it is seen from (22) that 


(27) |P|+1Q|+1S|=p-p6+6+0g=1-q™ <1. 


Thus, if (25) and (26) are satisfied, the largest of the sums of absolute values of the 
elements of the rows of matrix A is equal to 1, and therefore the difference scheme 
under consideration is stable. 

The inequalities (25) and (26) may be written 











nt < 1 
(Ax)? ~ 2 
nt (b + c)At | vast 1 | b+e 
£_ = 
2 a = o 5 ae o i+ bar * ; o 41) 20 
whenee 
2no 
2 \) > 
(28) (Ax) = er 
, A?’ o 
< - ae 

(29) ws min| |. 


In conclusion, choosing the intervals Ar and Af in accordance with (28) and (29) 
will insure the stabilitv of the difference scheme under consideration. 
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Finally, consider the problem of heat conduction in a cylinder subject to the 
boundary conditions stipulated in the second of Reid’s articles mentioned above. 
If, for the sake of convenience, we denote by z the radial distance ordinarily 
denoted by r, the mathematical formulation of the problem is as follows: 


oT oT 10aT 
a) ae +12) 
(31) T(x,0) = f(x) 
(32) ct te ome 
Ox 
(33) _K (37) = m(T(a,t) — v(1)] 
Ox z=a 
(34) Qmpo Co d » m, [T(a,t) — v(t)] — me v(t). 


ot 


The only essential difference between the above problem for the cylinder and the 
corresponding problem for the slab is the differential equation (30) which takes the 
place of (10) and the new condition (32). The difference analogs of (30) and (32) are 
| ee bes -_ = k | 7 — 2S vim + > grey 1 . Tete xt Tents] 
Ax 


At (Ax)? + mAx 2 








or 


in 1 i. 
(35) Tames =f ( = x) ae + (1 aa, oo + r(1 + a) Tutte 
m = 1,2,3,---M 
and 


(36) at = lau ° 


The difference equation (35) takes the place of the difference equation (17). In view 
of (35) and (36) it is readily seen that the previous developments for the slab may 
be carried out with the sole exception that in lieu of (24) the matrix A is now 


1 — 2r a" | 


1 1 
(1 _ i} 1 — 2r (1 ~ 1} 
1 


1 1 
roe a =a — 2 ——— 
(1 say 1 — 2r (1 + » 


% P Q S 


(37) A 


Il 
ae 
— 
| 

oI 
al 
— 
| 
bo 
bee 
as, 
os 
+ 
ol-— 
ee 

7 
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where the symbols P, Q, and S have the same significance as before. Since for the 
first M rows of A the largest of the sums of the absolute values of the elements is 
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equal to 1 (provided r < 4) we reach the conclusion that the criteria of stability 
are identical with those of the previous problem and are given in the inequalities 
(28) and (29). 

It may be briefly mentioned that the above analysis of stability may be readily 
extended to the case where the cylinder and the coaxial thin cylindrical shell are re- 
placed by a sphere and a concentric thin spherical shell. 
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New York, New York 


1. A. N. Lowan, “‘Heat conduction in a solid in contact with a well-stirred liquid,’ Phil. 
Mag., v. 17, 1934, p. 849-854. 

2. WauTER P. Ret, “‘Heat flow in a half space,” Quart. Appl. Math., v. 14, 1956, p. 206-208. 

3. Wa.LTER P. Rerp, “Heat flow in a cylinder,” Quart. Appl. Math., v. 16, 1958, p. 147-153. 

4. A. N. Lowan, The Operator Approach to Problems of Stability and Convergence of Solu- 


tions of Difference Equations and the Convergence of Various Iteration Procedures, Scripta Mathe- 
matica, New York, 1957, p. 79. 





On Numbers of the Form n‘ + 1 


By Daniel Shanks 


1. The Number of Primes. Let Q,(N) be the number of primes of the form 
n‘ + 1for1 < n S N. By a double sieve argument similar to that used for primes 
of the form n’ + a, [1], and for Gaussian twin primes, [2], one is led to the follow- 
ing conjecture: 








1 "dn 
(1) Q(N)~F 4 | log n 
where 
—1 2 —2 

(2) fil: G)+@)+G)), 

p=3 7 1 
the product being taken over all odd primes with () the Legendre symbol. Now 

s L (1) L2(1)L-2(1) ( 91 (239) 
= 1 — -&— 

* §°(2) wal, pP/\p-1 


where this product is taken over all primes of the form 8m + 1 and L,(s) and ¢,(s) 
are as defined in [1, p. 323]. We may therefore write 


(4) o* 328 a /2) wl. (3 a ;) (525). 


To evaluate this slowly convergent product we use the identity 


@ 1-4 CeIGE NO Cs ~ 
1+2 1+ 2 1+ 2 1+ 2 
which is valid for x < 4, and the identity 
2 8 2 
«) NOTED EAD _ L's (z + ‘) ; 
which is valid for s > 1. From tables of {.(s) and L.(s) we thus obtain 
(7) 8, = 2.67896 --- 


and therefore 














dn 
log n~ 





(8) Q(N) ~ Q(N) = 0.66974 [ 


It is interesting to compare this formula with that for the conjectured number 
{1] of primes of the form n’ + 1, 


N 
(9) P\(N) ~ P,\(N) = 0.68641 | :.. 
2 logn 
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TaBLe 1 
N | uw) | aww) Ovds 
100 | 18 | 19.5 | 0.924 
200 30 32.9 | 0.911 
300 | 44 45.1 0.976 
400 | 52 | 56.5 0.920 
500 63 67.5 0.934 
600 75 78.1 0.960 
7 | 80 88.4 0.905 
800 | 94 98 .6 0.954 
900 | 98 108.5 0.903 
1000 | 109 118.3 0.922 





The coefficients are nearly equal and have analogous formulae: 








4 —1 
“coapele Gg) 
0.68641 = - IT k se 
(10) w ; 
et G++) 
0.66974 = ; IT k = at uA 


2. A Table. A comparison of Q,(N) with the actual counts Q,(N) is handicapped 
by the very rapid increase in n‘ + 1. The 109th prime is already 984 095 744 257, 
nearly a trillion. A. Gloden [3] has completed the factorization of all n‘ + 1 up to 
n = 1000, following the work of Cunningham and others. He has kindly counted 
the primes for us, where 400 < n < 1000, and using his results we present Table 1. 
The deviations of Q,/Q, from unity are not unduly large considering the relatively 
small upper limit for N. For P;(N) and for the ordinary prime count x(N) we have 
similar deviations for N = 1000; x(1000)/l(1000) = 0.951 and P,(1000)/ 
P,(1000) = 0.924. 


3. Four Classes of Numbers. When we consider that Euler determined P,(N) 
up to N = 1500 over two hundred years ago [4], the present table of Q,(N) up to 
N = 1000 seems rather meager. The much greater difficulty of factoring the n* + 1 
numbers is fundamentally due to their much greater magnitude—but there are 
interesting technical differences also. The sieve method for n* + 1 used by Gloden, 
Cunningham, and others has three phases. 

A. Compile a list of primes of the form 8m + 1 

B. For each such prime solve the congruence 


x’ = —1 (mod p) 
“ <p : 
for its four roots. (Given one solution 2, the remaining three are congruent to —2, 
a, and —2;°.) 
C. With each x and each p divide out a factor of p for each n = x; + mp. 
Similarly determine those n‘ + 1 divisible by p’, p’, ete. 
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Now unfortunately there is much waste computation here. For instance, the 
hundred n‘* + 1 for n < 100 have 122 different prim-2s of the form 8m + 1 as factors. 
Yet all 295 of the 8m + 1 primes <100° must be examined in phases A and B, since 
a priori any such prime may be a factor of the n* + 1. And clearly this waste in- 
creases rapidly with N,—for N = 1000 we must examine all 19552 of the 8m + 1 
primes <1000° to factor out the (approximately) 1300 distinct actual prime 
factors. 

On the contrary, in the author’s sieve [5] for n’ + 1 there is no waste computa- 
tion and no phases A and B, either. The primes arise automatically in the sieve it- 
self, together with the corresponding solutions of the congruence, x” = —1 (mod p). 

This significant difference comes about as follcws. For every n, n* + 1 either 
has no new prime factor (n is “reducible”’) or it has precisely one new prime factor— 
and that to the first power (n is “‘irreducible”’). Therefore, if all prime factors cor- 
responding to smaller values of n have already been sieved out, each new prime 
stands exposed at the smallest n which satisfies n? = —1 (mod p). But for n* + 1 
we have not two but four classes of n; there are either 0, 1, 2, or 3 new prime factors 
in n‘ + 1. It is the occurrence of the “double” and “triple” irreducibles (i.e., 2 and 
3 new primes) which prevents the use of the automatic, n* + 1 type sieve for n‘ + 1. 
Already for n = 10 we have a double irreducible 


10° + 1 = 73-137, 


with the two new primes 73 and 137. 

Let R(N), Th(N), I2(N) and I;(N) be the number of “reducibles” (no new 
prime) and single, double, and triple irreducibles respectively which are < N. For 
example, 7,(120) = 92 and J,(120) = 28. Further, R(120) = J;(120) = 0, since 
neither reducibles nor triple irreducibles arise for nS 120. For larger n (from Glo- 
den’s tables) we find both reducibles 


29588 + 1 = 177-41-113-1249-16073-28513 
and triple irreducibles 
23762* + 1 = 637489-693569 -721057, 


but they are rare. 
The mean number of new primes is 





(11) (N) = BON) + =e) + 31(N) 


and in analogy with the situation for n’ + 1 the question arises whether »(N) has 
a limit for N — ©. For n? + 1, John Todd [5, p. 83] has conjectured »(N) > 
log 2 = 0.693. For n* + 1 and a modest N we have »(N) ~ 1.3. Analogy with 
Todd’s results concerning n® + 1 and log 2 would suggest a limit of log 4 for n‘ + 1, 
but there is no serious evidence in favor of this. 
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TECHNICAL NOTES AND SHORT PAPERS 


On Modular Computation 


By Henry B. Mann 


It will be assumed that the reader is familiar with the elementary theory of con- 
gruences [1]. 


Let m,, m2, --- , m, be s integers relatively prime in pairs and let M = mm, 
-++m,. Let 2%, %2, +--+ , 2 be an ordered set of s integers such that 0 < 2; < m;. 
There exists one and only one residue class x mod M such that z = x; (mod m,) 
and we therefore write z = (1 ,%2,-°-- ,%).Ifx= (m,--- ,%s),y = (Wi, °-* 5 Ys) 
then + y = (MH%,°°',%s + Ys), TY = (U1, -*- , Lee), Where the 7th co- 
ordinates must be reduced mod m;. The symbols (2, , --- , x,) are called modular 
numbers, but it should be kept in mind that (x, --- , z,) does not denote a num- 


ber but a residue class mod M. 

The purpose of this note is to describe a simple iterative procedure to determine 
the least non-negative residue mod M of a given residue class (2 , --- , 2). 

The iteration process described below gives the least non-negative residue in 
mixed radix representation. 

Notation. The moduli are denoted by m,, --- , m,. We put m = 1, 


M; = (0, --- , 0, 1,0, --- , 0); t#=1,---,8 

where 1 occurs in the ith position, 
. rT, = Mmm +--+ M4 ~=1,---,8, 
[x] denotes the least non-negative residue corresponding to the modular number 


z= (%,°°° ,2e)- 
The representation 


—1 
[elw = Qeaws, OS a;<m;, t=1,---,(8—-—1) 


is the mixed radix representation of [x]y . Because of the condition 0 < a; < m;, 
the mixed radix representation of a number a can be coded as a modular number 
La = (a, ey , Gs). 

Now proceed as follows: 

Find 


M; = Do aiwn;, x; = >. 2;;M; (mod M). 
3 


The jth columns of the matrices (a;;) and (2;;) are residues mod m;, so that 
the rows of these matrices may be regarded as modular numbers. 

In forming (2, --- , %,)(a:;) we proceed as in ordinary matrix multiplication 
but compute the inner product 2,a,; + --- + 2,a,; mod m;. Similarly, in forming 
(a; , -*- , @,)(x;;) we compute the inner product a,27,; + --- + a,2,;mod m;. Note 
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that if [zlyw = >> ajr;, then (a, --- , a,)(2;;) = 2, but z(a;;) gives in general 
not the radix representation of z, but of (2(a;;) )(24;)*. 
Given zt = (m, --~ , Z,), the iteration proceeds as follows: 


xz Az Aa a 
r= (t1,-* 5X) - © = 2(a;;) 
zr" 2A a‘ (2:5) 2s a x (x a x) (a;;) a” ae a” + Aa” 
ry” =a a‘” (x;;) os zg” (x a x) (a;;) a‘**» a a” + Aa‘ 


The process ends when Ar = 0. 

The first two coordinates of x” will coincide with the first two coordinates of z. 
Thereafter the first « + 1 coordinates of x“ will coincide with the first (a + 1) 
coordinates of x so that x“ = z and hence a®~” gives the mixed radix representa- 
tion of [2] , the coordinates of a” being the digits of this representation. In many 
cases, however, less than (s — 1) steps will suffice. 


Example 1 (ascending order): m, = 2, m2. = 3, mz = 5, m = 7. 


a a a a 
PM, 212 ae 22 2 
wn OEP eg wig Hag 
00 0 4 000 2 
z Az Aa a 
1, 2, 2,5 1,2, 1,0 
m= 1, 2,1, 4 6, &.4,..1 0, 0, 1, 1 1 2,:3,4 
% = 1, 2, 2,5 
Hence [zla0 = 1 + 2.2 + 2.6 + 1.30 = 47. 
Example 2 (descending order): m, = 7, m2 = 5, ms = 3, m = 2. 
re Su9*t Be “ore 
1 Ge Te es 1k a are 
FE Qe 86 es Fp Bs 
o- 0° 6" 4 000 1 
Zz Az Aa a 
4, 3, 1,0 2.2.3 
ra ee 0, 0, 0, 1 0, 0, 0, 1 4, 2, 2,0 
4, 3,.1,.0 


? ? 


Hence [4, 3, 1, Oleo = 4 + 2.7 + 2.35 = 88. 


> 


Remarks: (1) Arranging the modules in ascending order may have an advantage 
because the matrices (a;;) and (x;;) are half diagonal and if 0 < y; < m;thenalso0 < 
yj; < ™j4, so that y; never has to be converted. However, any arrangement will 
work. 

(2) No conversion to decimals is necessary to decide if [z’]y > [2°]. If 
(ele = do; 0;x;, = 2s a;”x;, then [x”|y > [2] if and only if 

(2) “a (2) () (2 
a,° > 4, OF as i one | @5-448 = As-ju1, aS”, > Gs-;. 


* Caution: The matrix multiplication defined above is not associative. 
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(3) All calculations are mod m; so that the digits of the mixed radix representa- 
tion of [z]4 can be obtained using only calculations mod m; . 
(4) The value x — x can be obtained from x — 2 = x — x + 2°” — 
z” so that the modular number z need not be remembered during the whole 
process. 

(5) The matrices (a;;) and (2z;;) are computed preliminary to the iteration 
procedure and are not part of it. 


a) 


Ohio State University 
Columbus 10, Ohio 


1. B. M. Stewart, Theory of Numbers, Macmillan, 1952, p. 111-113 and p. 130. 


Generation of Permutations by Transposition 
By Mark B. Wells 


1. Introduction. As discussed by Tompkins [1], many problems require the 
generation of all ~! permutations of n marks (henceforth called arrangements). 
This note presents a generation scheme whereby each step consists of merely trans- 
posing two of the marks. The bookkeeping is quite simple, thus this scheme is some- 
what faster than either the usual dictionary order method or the Tompkins-Paige 
method [1]. Also, the important property of leaving the (j + 1)st position alone 
until all 7! arrangements of the marks in the first 7 positions have been generated 
is preserved. 


2. Notation. An arrangement of n marks will be given by an n-tuple, (m,, mz, - - - 
m,). A permutation, that is, an operation of permuting an arrangement of marks, 
will be given in cyclic form, with P’s modified by subscripts as entries. The subscripts 
indicate the position of the marks to be moved in the n-tuple on which the permuta- 
tion is operating. For example, if a = (1, 2, 5, 4, 3) is an arrangement of five marks 
and p = (P;P;P2)(P.Ps) is a permutation, then p(a) = (2, 5, 1, 3, 4). 

The bookkeeping for this generation scheme is handled, as in most schemes of 
this type, by an ordered set of indices ¢, , k = 2, 3, --- , n, where each & assumes 
the values 1 through k and indicates the progress of the subgeneration of the ar- 
rangements of marks in positions 1 to k. (This is essentially the “signature” dis- 
cussed in [1].) Thus there are n! sets of values for the ¢,’s, one set for each arrange- 
ment of the n marks. The set ¢, = 1 for all k corresponds to the initial arrangement, 
and successive sets are formed in dictionary order (assuming increasing significance 


with increasing subscript). An index k’ gives at each step the smallest subscript k 
for which t, ¥ k. 


3. The Generation Rules. The transposition required at each step depends on 
the current value of the index k’ and on the corresponding value of ty: 4:(ta4: is as- 
sumed = i). The rules are: 

I. If k’ is even, then interchange the marks in positions k’ and k’ — 1. 
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II. a. If k’ is odd and t-4, < 2 then interchange the marks in positions k’ and 


k’ — 1. 
b. If k’ is odd and Z < t+4;, < k’, then interchange the marks in positions k’ and 
KW — tera + 1. 


ce. If k’ is odd and &-4, 2 k, then interchange the marks in positions k’ and 1. 

Before proving that these rules yield all n! arrangements in n! — 1 applications 
(starting with a given arrangement), let us illustrate their application for n = 5. 
The rules apply at step s to yield the arrangement given at step s + 1. 

















Step te t te t & k’ » oe vr. 
1 OS om ory 2 (1, 2, 3, 4, 5) 
2 si. 5 3 (2, 1, 3, 4, 5) 
3 (2 PRs 2 (2, 3, 1, 4, 5) 
4 22111 3 (3, 2, 1, 4, 5) 
5 bot bs od 2 (3, 1, 2, 4, 5) 
6 $82 Ava 4 (1, 3, 2, 4, 5) 
7 Br ER 2 (1, 3, 4, 2, 5) 
12 33's 4 4 (i, 4, 3, 2, 5) 
13 P'S 2 (1, 4, 2, 3, 5) 
18 $18s88.2 4 (2, 4, i, 3, 5) 
19 Prnn 2 (2, 4, 3, 1, 5) 
24 fh oe ee 5 (3, 4, 2, 1, 5) 
25 ek 2 (3, 4, 2, 5, 1) 
48 23421 5 (2, 5, 4, 3, 1) 
49 e 4ute Bs 2 (2, 5, 4, 1, 3) 
72 ae a 5 (4, 1, 5, 2, 3) 
73 1 1 4 1 2 (4, 1, 5, 3, 2) 
6 | 23441 5 (5, 3, i, 4, 2) 
97 | 111651 2 (5, 3, 1, 2, 4) 
120 | $3451 | 6 (i, 2, 3, 5, 4) 





A close inspection of the above example will reveal the mechanism at work. 
Following a transposition (P;P;,) with i < k all (k — 1)! arrangements involving 
change only in positions 1 through k — 1 are generated before P, appears again. 
During a complete subgeneration of the k! arrangements of the k leftmost positions, 
the transposition (P;P,), for some particular i < k, occurs k — 1 times, each time 
k’ = k. The particular value of i will be k — 1, k — t4; + 1, or 1, according to 
the rule in force. To insure that no duplicate arrangements appear, the mark 
initially (at the time the subgeneration begins) in position k and the marks suc- 
cessively (each time (P;P;) is performed) in position 7 must all be distinct. This is 
accomplished in two ways according as k is even or odd. For an example with k = 4, 
compare the marks in position 4 at step 1, and in position 3 at steps 6, 12, and 18. 
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Lemma 1. Let a = (P,P yoy Px) with i, , t2, = ste — k- 1 bea cycle and 
let 7 < k. Then of(P;P,)al** = (P;Px). 

Proof. This is verified by direct permutation multiplication. 

The significance of this lemma is the following. Let k be odd and consider any 
subgeneration of k! arrangements of the marks in the k leftmost positions. During 
this subgeneration k’ will be equal to k k — 1 times, and we will have k — 1 identi- 
cal applications of rule II interspersed with k identical permutations of the first 
k — 1 positions. If, as Lemma 2 will show, this permutation is a cycle, then Lemma 1 
says the effect of the entire subgeneration was as a single application of rule IT on 
the initial arrangement. 

Lemma 2. For k even, (P:Ps-1)(PsaPx) ( []i=i (PsPea)(PiaPr) (PrP) = 
pm, @ single cycle, where pp = (PiP2), px = (PiPsP2P3) and in general, 
ed (PiPiPx-ol PrsPis aes P3} Peal Pe+Pi-s dyn P3}). 

Proof. Again, direct multiplication gives verification. 

Thus for k even the effect of complete subgeneration is to permute the k marks 
by a cycle. For examples of the effects given by these two lemmas, compare the ar- 
rangements at steps 1 and 24 and at steps 25 and 48 (k = 4) and at steps 1 and 120 
(k = 5). 

Consider now any such subgeneration beginning with the arrangement, say 
(m , M2, +++ , 7m, °*- m,). With k even, each of the k — 1 applications of rule I 
finds a new mark to put in the kth position, since during this subgeneration ¢, is 
assuming the values 1, 2, --- k, and so by the special construction of rule II (and 
Lemma 1), position k — 1 contains successively, at the time of application of rule 
I, mMy—2 , M1, Me—zs , Mes, -** , Mm. With */ odd, each of the k — 1 applications of 
rule II finds the k — 1 leftmost marks permuted by a (k — 1)-cycle (by Lemma 2), 
and hence also finds a new mark to put in the kth position. A simple induction now 
shows that such subgenerations yield k! distinct arrangements. We have proved 
the following: 

THEorREM. The generation scheme as given above yields all n! arrangements of n 
marks in exactly n! — 1 steps (starting from a given arrangement). 


4. Remarks. As in most other generation schemes the property of changing the 
jth mark only when all arrangements of the previous 7 — 1 marks have been gener- 
ated allows significant time-savings in some problems. If following a transposition 
(P;Pi41) witht < k + 1 the problem decides it does not need to use the k! arrange- 
ments formed by permuting the present k leftmost marks, then this subgeneration 
may be skipped by applying Lemma 1 or Lemma 2 according as k is odd or even. 
This immediately prepares the arrangement for the next application of (P;P.4:). 
The permutation of Lemma 2 is not a transposition, but is quite easy to code into 
the scheme. An interesting question is whether or not an equally simple generation 
by transposition scheme exists in which block skipping is also always done by 
transposition. 

With the assumption that the marks being permuted are in most problems indices 
used for address modification, and thus should occupy the address portion of a 
computer word, a time comparison [2] between this scheme and the Tompkins-Paige 
method was made on Maniac II. With nine marks the transposition scheme gener- 
ates arrangements about twenty per cent faster. In addition, the transposition 
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scheme is advantageous for problems in which minimum mixing of the marks at 
each step is important. 
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Chebyshev Approximations to the Gamma 
Function 


By Helmut Werner and Robert Collinge 


In this note several Chebyshev approximations are given for the function 
= I(x + 2) for x in the 0 S x S 1.0 range. The approximations were obtained 
from a table of !'(z + 2), employing well-known methods as described in numerous 
papers; see for instance [1] and the literature quoted there. The table of T(z + 2) 
was calculated from the asymptotic expansion of log ['(z) as given in [2] to provide 
data accurate to at least 10". Compare also [3]. 
The asymptotic expansion of In ['(z) is given by 


In T(z) = (2 — 4) nz — 24+ In V2 + ®(z) 
where 
@(z) = > $ 1m A. + R,(z), 
= 2r(2r — cork = 


and B, is the rth Bernoulli number. 

It can be shown [2] that for z > 0 the value of @(z) always lies between the sum 
of n terms and the sum of (n + 1) terms of the series, for all values of n. In ter- 
minating this series with the nth term the error R,,(z) will be less than 


0" ag cee eee 
2(n + 1)(Qn +1) 2° 
By truncating @(z) at the 10th term it is easily shown that for values of z 2 13, 


the error in the expansion is less than 5.5 X 10°”. We therefore replace 6(z) by 
> 2: A,/2** and calculate In ['(z) for values of z in the range 13 S z < 14. 


Received July 25, 1960. 





TABLE 1 





10 





= 
a 
2 


f 
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Table of Coefficients 
7 & 

0.25 X 107 0.16 X 10° 
0.99999 99758 0.99999 99998 452 
0.42278 74605 0.42278 43662 730 
0.41177 41955 0.41183 92935 920 
0.08211 17404 0.08159 03449 474 
0.07211 01567 0.07416 00915 535 
0.00445 11400 0.00007 55964 181 
0.00515 89951 0.01033 20685 065 
0.00160 63118 | —0.00157 80074 635 

0.00079 62464 760 


0.99999 99999 9269 
0.42278 43369 6202 
0.41184 02517 9616 
0.08157 82187 8492 
0.07423 79076 0629 
—0.00021 09074 6731 
0.01097 36958 4174 
—0.00246 67479 8054 
0.00153 97681 0472 








—0.00034 42342 045 
0.00006 77105 711 





6 
7 





13 


15 





fe! » 


0.96 X 10° 


0.97 X 10-* 





= 


WOWONMAH or WN © 


0.99999 99999 99990 44 
0.42278 43351 02334 79 
0.41184 03301 66781 29 
0.08157 69261 24155 46 
0.07424 89154 19444 74 
— 0.00026 61865 94953 06 
0.01114 97143 35778 93 
—0.00283 64625 20372 82 
0.00206 10918 50225 54 
—0.00083 75646 85135 17 
0.00037 53650 52263 07 
—9.00012 14173 48706 32 
0.00002 79832 88993 83 
—0.00000 30301 90810 28 


0.99999 99999 99999 9032 
0.42278 43350 98518 1178 
0.41184 03304 21981 4831 
0.08157 69194 01388 6786 
0.07424 90079 43401 2692 
—0.00026 69510 28755 5266 
0.01115 38196 71906 6992 
—0.00285 15012 43034 6494 
0.00209 97590 35077 0629 
—0.00090 83465 57420 0521 
0.00046 77678 11496 4956 
—0.00020 64476 31915 9326 
0.00008 15530 49806 6373 
—0.00002 48410 05384 8712 
0.00000 51063 59207 2582 
—0.00000 05113 26272 6698 





17 


18 





a 
= 


0.10 X 10-7 


0.10 X 1078 





OCOWONAUrPWNHH © 





0.99999 99999 99999 99901 2 
0.42278 43350 98467 79580 6 
0.41184 03304 26367 20638 1 
0.08157 69192 50260 90508 9 
0.07424 90106 80090 41696 9 
—0.00026 69810 33348 38176 8 
0.01115 40360 24034 39169 2 
—0.00285 25821 44619 65607 6 
0.00210 36287 02459 83329 2 
—0.00091 84843 69099 08014 2 
0.00048 74227 94476 75810 4 
—0.00023 47204 01891 94985 9 
0.00011 15339 51966 59947 0 
—0.00004 78747 98383 44672 4 
0.00001 75102 72717 90508 0 
—0.00000 49203 75090 42313 2 
0.00000 09199 15640 71621 4 
—0.00000 00839 94049 59039 7 





0.99999 99999 99999 99990 
0.42278 43350 98467 21319 
0.41184 03304 26430 62304 
0.08157 69192 47528 84581 
0.07424 90107 42094 91715 
—0.00026 69818 88740 38315 
9.01115 40438 29069 91793 
—0.00285 26318 64702 11862 
0.00210 38579 20672 20524 
—0.00091 92675 95039 95026 
0.00048 94361 06998 14458 
—0.00023 86428 33752 63647 
0.00011 73283 10224 09396 
—0.00005 43183 86280 13508 
0.00002 28140 41153 66022 
—0.00000 80523 43363 48309 
0.00000 21741 77495 45532 
—0.00000 03889 70057 38769 
0.00000 00339 81801 01810 


02 
64 
23 
87 
38 
07 
28 
89 
09 
11 
34 
10 
51 
99 
75 
46 
64 
55 
43 
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For the convenience of the reader the A; coefficients are quoted below, to 25 
significant figures. 


A, = 0.08333 33333 33333 33333 33333 3 

Ag = —0.00277 77777 77777 77777 77777 78 
A; = 0.00079 36507 93650 79365 07936 508 
A, = —0.00059 52380 95238 09523 80952 381 
As = 0.00084 17508 41750 84175 08417 508 


Ag = —0.00191 75269 17526 91752 69175 27 
A; = 0.00641 02564 10256 41025 64102 56 
As = —0.02955 06535 94771 24183 00653 6 


Ay = 0.17964 43723 68830 57316 49385 
Aw = —1.39243 22169 05901 11642 7432 
In+/2x = 0.91893 85332 04672 74178 03297 


A triple precision logarithm routine was used to evaluate In z, and then an ex- 
ponential routine to calculate T(z) = e'"". Each of these routines produces re- 
sults accurate to at least 24 significant digits. 

After obtaining a table of ['(z) for z in the range 13 S z S 14, we made use of 
the recursion formula T'(z + 1) = zI'(z) in order to obtain a table of T(z + 2) for 
zx in the range 0 S x S 1.0. 

From the tests made on the results obtained, the values of [T(z + 2) were 
shown to be accurate to at least 21 significant figures. 

Several Chebyshev approximations have been calculated to provide varying de- 
grees of accuracy. Let 


r(2 +2) = Diaz" + «(z) 
v=) 


and 


6. = max | €n() |. 
Oszs1 


Table 1 gives the coefficients a,“” for n = 7, 8, 10, 13, 15, 17, 18 together with the 
corresponding «42 . 


Professional Services Department 
Burroughs Corporation 

460 Sierra Madre Villa 

Pasadena, California 


1. E. L. Streret, “Numerical methods of Tchebycheff approximation,’’ On Numerical 
Approximation by R. E. Langer, The University of Wisconsin Press, Madison, 1959, p. 217-232. 

2. E. T. Wuirraxer & G. N. Watson, A Course of Modern Analysis, Cambridge, 1958, p. 
251-253. 

3. M. E. SHerry & S. Futpa, “‘Calculation of gamma functions to high accuracy,’’ MTAC, 
v. 13, 1959, p. 314-315. 

4. H. 8. Unuer, “Log =z and other basic constants,’’ Proc. Nat. Acad. Sci. U. 8. A., v. 2A, 
1938, p. 23-30. . 

5. H. S. Unter, “‘The coefficients of Stirling’s series for log I'(x),’’ Proc. Nat. Acad. Sci. 
U.S. A., v. 28, 1942, p. 59-62. 
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A Partition Test for Pseudo-Random Numbers 
By J. C. Butcher. 


A frequently used test on random digit sequences consists in forming from them 
sets of n integers; 1, 2, --- , vn say, which, on the hypothesis of randomness, are 
independent and have equal probabilities, 1/z, of being each of the integers 1, 2, - - -, 
x (x usually being chosen as a power of 2). From any equalities that may exist be- 
tween 1, %, °° , ¥ & partition of n is defined and the actual test is on the fre- 
quency of occurrence of the different partitions of n. 

For example, a popular form of the test known as the “‘poker test”’ distinguishes 
the different partitions of n = 5 known by the descriptive names “all different,”’ 
“one pair,” “two pairs,” “three of a kind,” “full house’’ (three of one kind, two of 
another), “four of a kind” and “five of a kind.” 

It is obvious that much computing time must be absorbed in distinguishing be- 
tween the various possibilities, and that if n is much greater than 5, the number of 
possible partitions becomes very large. For this reason it may be advantageous to 
group the partitions together in some way that lowers the number of cases and also 
simplifies the programming necessary to distinguish these cases. 

A convenient way of doing this will now be described. It has the added advan- 
tage that the calculation of the expected probabilities is extremely simple. Let us 
distinguish n classes of partitions C, , C., --- , C, where C, includes all partitions 
defined from exactly r different integers occurring amongst », , v2, --- , ¥ . In the 
example n = 5, C, would include the single partition “five of a kind,”’ C2 would in- 
clude the two partitions “full house” and “four of a kind,” C; would include “three 
of a kind” and “two pairs,” while C, and C; would each contain the single partitions 
“one pair” and “all different,’ respectively. 

If the n integers 1; , v2, --- , vy, are generated one after the other, and if we define 
Si = 1, 2, --- , n) as the set of integers 1, 2, --- , z, other than those identical 
with v; for some 7 < 7, then the value of r that characterizes the class C, is given 
by 


r= re €i, 
i=l 
where ¢; is equal to 1 or 0 according as »; is or is not a member of S; . 

When z is no more than the computer word length, then S; is conveniently rep- 
resented by the binary digits of a single word, and this word is quickly computed 
from the word representing S;_, if, for example, a logical conjunction instruction is 
available. 

To find the expected frequency of occurrence of the class C, for r = 1,2, --- ,n 
we must find the probability p, that the class C, will arise in a single trial. It is easy 
to see that 
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(z—1), 


gn 2% 


(x — 1)(@ — 2)----(2 —n + 1) 
a, 


Pn es gz 





, 


where a; , @2, --- , @, are dependent on n, but not on z. Since p, + pe + --- + 
Pn is identically equal to unity, we find 
(1) x" = ar + ae(e — 1) + --- + aar(x — 1) --- (@ —n +1) 


for all values of x. The coefficients a; , a2, --- , a, are thus identical with Stirling’s 
numbers of the second kind of order n [1], and they can be evaluated in turn by 
substituting = 1, 2, --- , in (1) and using as a check the value a, = 1. 


Department of Applied Mathematics 
University of Sydney 


1. C. Jornpan, The Calculus of Finite Differences, Chelsea Publishing Co., New York, 1947, 
p. 170. 





REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


34(D, E, L, S, X].—Nextson A. Locan, General Research in Diffraction Theory, 
Volumes I and II, Missiles and Space Division, Lockheed Aircraft Corp., 
Sunnyvale, California. Reports LMSD-288087 and LMSD-288088, December 
1959, xxiii + 345 p. and xviii + 268 p., 28 cm. Deposited in UMT File. 


Volume I consists of a study of the theory and applications of a class of inte- 
grals defined by 


An’(—) = Yan [Ai(— ay) Meme, 


Bn" () = > B."Ai( —B,) et 


where a,, 8, denote the roots defined by Ai(—a,) = 0, Ai’(—8,) = 0, and 
Ai'(—a,), Ai(—8,) are the turning values of the Airy integral. This representation 
is useful only for ¢ > 0. Alternative representations useful for § — 0 are developed 
for the case Ao” and By". For m = 1 the functions are entire functions of £, and 
tables are given for the coefficients of the Taylor series of A," and B,”". These coeffi- 
cients are evaluated by using the Euler summation scheme to sum the divergent 
series obtained by setting = 0, m = 1 in the above representations. When m = 2 
it is necessary to extract some terms which are singular at = 0. The remaining 
parts of A,” and B,” are shown to be entire functions. Tables for the coefficients 
in the Taylor series for these non-singular parts are found by using the Euler-Mac- 
laurin summation formula to sum the divergent series which are obtained by setting 
& = 0, m = 2 in the above representations. For  — — ©, asymptotic expansions 
are obtained for the cases m = 1 and m = 2. Tables are given for the coefficients 
in these expansions. 

Volume II consists of a set of 26 tables and 17 figures that provide a supplement 
to the theoretical analysis contained in Volume I. Tables A, B, and C contain 
special tables of exponential and trigonometric functions which will facilitate com- 
putation with residue series and asymptotic expansions of the diffraction integrals. 
The functions tabulated in the remaining tables can be used to study diffraction 
effects when (a) source and receiver are on the surface, (b) source (or receiver) is 
on the surface and the receiver (or source) is at a great distance, and (c) both source 
and receiver are at a great distance. 

AUTHOR’s SUMMARY 


35[F].—_R. Kortum & G. McNrgt, A Table of Quadratic Residues for all Primes 
less than 2350, LMSD 703111, October 1960, Lockheed Missiles and Space 
Division, Sunnyvale, California, iii + 3 + 378 unnumbered pages, 28 cm. 


This large report, bound with a plastic spiral, lists all 187,255 quadratic residues 
of the 347 primes from 3 to 2347. The tables were computed on an IBM 7090 in 
about ten minutes. Presumably most of this time was spent in binary-decimal con- 
version and in writing on tape. The original printing was done on a high-speed, wire 
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matrix printer and is readable, but certainly not elegant. In compensation, the 

tables are very easy to use, since the spiral binding allows the pages to lie flat. 
The tables also give N(p), the number of positive non-residues <}$p. In the 

introduction it is pointed out that for primes of the form 4m + 3 we have 


[3(p — 1)]! =(—1)*™ (mod p). 
It is also indicated that for all such primes (but we must add >3) the class number, 


h(—p), is given by 
p—l 
M—p)=->> (2)o. 
P a=1 


The much more easily computed formula [1], 


— i — 4N(p) 
h ~ ) = P ? 
(—?) = "T= FeP) 
is not mentioned. The introduction also states that it can be “found” in the table 
that N(p) = m for all primes of the form 4m + 1. But surely one does not need the 
table to be convinced of this simple theorem. The quantity which is really useful 
for those primes is 9>su (a/p), and not the redundant N(p). 





D. 5. 


1. E. Lanpau, Aus der elementaren Zahlentheorie, Chelsea Publishing Co., New York, 1946, 
p. 128. 


36[G X].—V. N. Fappeeva, Computational Methods of Linear Algebra, Translated 
by Curtis D. Benster, Dover Publications, Inc., New York, 1959, x + 252 p., 
21 em. Price $1.95. 


The first chapter of this book forms a clear and well-written introduction to the 
elementary parts of linear algebra. The second chapter deals with numerical meth- 
ods for the solution of systems of linear equations and the inversion of matrices, 
and the third with methods for computing characteristic roots and vectors of a 
matrix. Most of the important material in these domains is to be found here, and 
many numerical examples which illustrate the algorithms and point out their merits 
and deficiencies are given. 

The discussion is directed principally to the hand computer, and machine 
computation in the modern sense is hardly present, but the book must be regarded 
as a valuable guide for the worker in the general area of linear computation. 


Morris NEWMAN 
National Bureau of Standards 
Washington 25, D. C. 





37[G, X].—M. Mipuat J. Gazaus, Les Structures de Commutation a m Valeurs et 
Les Calculatrices Numériques, Collection de Logique Mathématique, Série A, 
Monographies Réunies par Mme. P. Fevrier, Gauthier-Villars, Paris, 1959, 78 p., 
24 cm. Price 16 NF. 


The theme of this pamphlet is sets of operations which are complete in the sense 
that “conjunction” and “negation” (or “exclusive or,” “conjunction” and “1,” or 
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the “Sheffer stroke’) are complete. By an operation we mean a (single-valued) 
function whose domain is S” for some set S and positive integer n. A set F of 
operations on S is complete if any operation f on S (of any number of argu- 
ments) can be constructed from F by composition (substitution) and identification 
of variables. 

The first three chapters, which are introductory, include, among other things, a 
discussion of how constructing f from F corresponds to constructing a network “‘re- 
alizing”’ f from elements (primitive networks) which realize the elements of /. Chap- 
ter IV is preparatory to Chapter V, where the first significant theorem appears. This 
is to the effect that sum modulo p (p a prime), product modulo p, and the constant 
functions are complete. Alternatively, every n-ary operation on 0, 1,2, ---,p— 1 
is representable by a polynomial in n variables over the field of integers, modulo 
p. The author fails to note, however, that for any finite field, any operation on the p” 
field elements is representable by a polynomial over the field. As a matter of fact, 
essentially the same argument the author gives for p elements is applicable to the 
more general situation. 

Chapter VI deals with a theorem of Webb to the effect that the binary operation 
W defined over 0,1, --- ,m— 1by W(a, y) = Oifz ¥ yand by W(2,y) = x +1, 
mod m, if x = y is complete. The author gives a formulation of the theorem which 
does not make use of the additive structure on the set, and gives a proof of it. 

The last chapter (VII) generalizes a theorem of E. L. Post to the effect that if R 
is a permutation of EZ, the integers mod m, then the pair of operations ®,z , Pe is 
complete, where R(i)@2R(j) = R (min i, j) and Pe(R(i)) = R(i + 1), for all 
jek. 

The following misprints were discovered: 

p. 39 line 10, E = (0,1 ---,” — 1) should read E = (0,1, ---,p— 1); 
p. 40 (63), read A, for Arq; 
p. 40 line 6 from bottom, read M,, for Mrt; 
p. 51 line 2 from bottom, read p” terms for p terms; 
p. 55, 56; each recurrence of bg should read b, ; 
p. 59 last line, a = dab should read qg = dab. 
Catvin C. ExGor 


IBM Research Center 
Yorktown Heights, New York 


38[I).—A. O. Ge.ronp, Differenzenrechnung, Deutscher Verlag, Berlin, 1958, viii 
+ 336 p., 23 em. Price DM 40. 


This is a translation of the Russian edition (1952) which is a revised and ex- 
tended version of an earlier book (1936). It is mainly concerned with problems in 
the complex domain, and some material, traditional in courses on Finite Differences, 
is omitted. The techniques used are those of classical analysis. There are occasional 
sets of problems, and some very interesting worked examples. 

The book is divided into three large chapters (1, 2, 5) and two smaller ones 
(3, 4). Chapter One deals with the problem of interpolation (‘“‘construct an (ap- 
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proximate) expression for a function, given its values at a discrete set of points”’), 
and includes an account of the Chebyshev theory. 

Chapter Two deals with the Newton Series, first for equally spaced nodes, then 
for more general cases. The chapter concludes with applications of interpolation- 
theoretic methods to number-theoretic problems, in particular, to a proof of the 
theorem that a and 8 = e* cannot both be algebraic, except for a = 0. 

The early part of Chapter Five is concerned with conventional material, includ- 
ing, for example, Ostrowski’s proof of Hélder’s result that T'(z) does not satisfy an 
algebraic differential equation; the latter part is concerned with work of the author 
(1951) on linear differential equations of infinite order, with constant coefficients. 

Chapter Three is concerned with earlier (1937) research of the author on the 
construction of (entire) functions given their values at a series of points a, ,a,— ~, 
and with related problems, e.g., the uniqueness of such functions. 

Chapter Four contains standard material on the Summation Problem and the 
theory of Bernoulli numbers and polynomials; it includes, e.g., a proper account of 
the Euler-MacLaurin Summation Formula. 

The book is clearly and precisely written. It can be recommended as an excel- 
lent source for many of the basic theorems in numerical analysis, and is a very 
suitable complement to such books as Natanson [1], which is largely concerned 
with the real domain. 

Joun Topp 


California Institute of Technology 
Pasadena, California 


1. I. P. Natanson, Konstrucktive Funktionentheorie, Akad. Verlag, Berlin, 1955 [MTAC 
Review 9, v. 13, 1959, p. 64-67.] 


39[I, X].—J. Kuntrzmann, Méthodes Numériques, Dunod, Paris, 1959, xvi + 252 
p., 25 em. Price NF 36.00. 


The author (who is a professor of applied mathematics of the Faculty of Sciences 
at Grenoble) admits his concern over the lack of a suitable textbook in numerical 
mathematics written in French. Rather than translate a foreign (to him) work, he 
decides to write a new book. 

For various reasons he decides to limit his book almost exclusively to interpola- 
tion. The usual interpolation formulas (Newton-Gregory, Stirling, Gauss, Bessel, 
Everett, and Lagrange) are included for equally-spaced abscissas and also for divided 
differences as appropriate. 

For the most part, approximation by the standard sets of polynomials (Legen- 
dre, Chebyshev, etc.) is avoided, but Bernoulli polynomials and Bernoulli numbers 
are discussed. 

More general formulas for which the given data might be either values of the 
function or values of certain derivatives are discussed. Numerical integration is 
avoided, but interpolation for functions ‘of two or more variables as well as of a 
complex variable is included. The last two chapters deal with the theory of inter- 
polation for linear sums of special functions (exponentials, trigonometric sums, etc.). 

Since the book was written to fulfill a need in France, and since there is no cor- 
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responding need in the United States, the reviewer feels that the book will have 
limited appeal to American numerical analysts. 


Rosert E. GrREENwoop 


The University of Texas 
Austin, Texas 


40[K].—W.S. Connor & Marvin ZELEN, “Fractional factorial experiment designs 
for factors at three levels,” Nat. Bur. Standards Appl. Math. Ser., No. 54, May 
1959, v + 37 p., 26 cm. For sale by the Superintendent of Documents, U. S. 
Government Printing Office, Washington 25, D. C. Price $.30. 


This is a sequel to NBS Applied Mathematics Series, No. 48 [1] which contains 
plans for fractional factorial designs for factors at two levels. The present compila- 
tion lists fractional factorial designs for factors at three levels as follows: for 1/3 
replications, 2 for 4 factors and 3 each for 5, 6, and 7 factors; for 1/9 replications, 
3 each for 6, 7, and 8 factors; for 1/27 replications, 3 each for 7, 8, and 9 factors; for 
1/81 replications, 3 each for 8 and 9 factors; and for 1/243 replications, 3 each for 
9 and 10 factors. For the same replication and number of factors the designs differ 
by the size of the blocks into which the treatment combinations are arranged. 
No main effects are confounded with other main effects or with two-factor in- 
teractions. Measurable two-factor interactions when the design is used as com- 
pletely randomized or when treatments are grouped into blocks are listed. In addi- 
tion, interactions confounded with blocks are given. Text material discusses the 
plan of the designs, loss of information, and the analysis of this type of designs. 

C. C. Craic 


The University of Michigan 
Ann Arbor, Michigan 


1. NBS Appurep Martuematics Series, No. 48, Fractional Factorial Experiment Designs 
for Factors at Two Levels, U. S. Gov. Printing Office, Washington, D. C., 1957 (MTAC Review 
7, v. 12, 1958, p. 66). 





41[K].—Epwin L. Crow & Rosert S. Garpner, “Confidence limits for the ex- 
pectation of a Poisson variable,’ Biometrika, v. 46, 1959, p. 441-453. 


For any m, the inequalities >°%2., pi(m) = 1 — ¢ 2pm) <1 — «¢ 
te" pm) <1 — efor all a, where p;(m) = e "m‘/i!, define c,(m) and c2(m) 
uniquely. Define mz(c) to be the smallest m for which c2(m) = c, and my(c) to be 
the greatest m for which c;(m) = ec. 

Table 1, p. 448-453, gives mz, and my to 2D for e = .2, .1, .05, .01, .001, and 
c = 0(1)300. The table was computed to 3D on an unspecified electronic computer; 
when the computed third place was a 5, the 5 was retained in the printed table. 

Table 2, p. 444, compares the present confidence limits with the system 6, of 
Pearson & Hartley [1] and the system 6; of Sterne [2]; table 3, p. 445, compares them 
with the approximate formulas of Hald [3]; table 4, p. 446, compares them with the 
mean randomized confidence intervals of Stevens [4]. 

Reprints may be purchased from the Biometrika Office, University College, 
London, W.C. 1, under the title “Tables of confidence limits for the expectation of a 
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Poisson variable.”’ Price: Two Shillings and Sixpence. Order New Statistical Tables, 
No. 28. 


J. ArtHur GREENWOOD 


Princeton University 
Princeton, New Jersey 


1. E. 8. Pearson & H. O. Hartiey, Biometrika Tables for Statisticians, Cambridge Univ. 
Press, 1954, p. 204-205. [MTAC, v. 9, 1955, p. 205-211]. 

2. T.E. Srerne, “Some remarks on confidence or fiducial limits,’’ Biometrika, v. 41, 1954, 
p. 275-278. [MTAC, v. 9, 1955, p. 216]. 


3. ANDERs Hap, Statistical Theory with Engineering Applications, John Wiley & Sons, 
prey | New Mer 1952, p. 718. 
W. . STEVENS, PsPiducial limits of the parameter of a discontinuous distribution,’ 
bicneehthe v. 37, 1950, p. 117-139. 


—F. G. Foster & D. H. Rees, Tables of the Upper Percentage Points of the 
Generalized Beta Distribution, New Statistical Tables Series No. 26. Reprinted 
from Biometrika, vols. 44 & 45. Obtainable from the Biometrika Office, Uni- 
versity College, London, W.C. 1, 30 p., 27 cm. Price five shillings. 


In sampling from a k-variate normal population, let A and B be independent 
estimates, based on », and » degrees of freedom, respectively, of the population 
variance-covariance matrix. Let 6, S --- < & be the roots of the determinantal 
equation | 64,A + (@ — 1)»B| = 0. Then the distribution of @, is given by 


z 6% Py A 
I,(k; p,q) = K | ao, [ d0,4°-- [ dé, II 6? "(1 — &)*" [] (6, — 4) 
0 0 0 t=1 i>i 


where p = }(| m — k| + 1), q = 4(m — & + 1). I,(1; p, g) is simply the incom- 
plete-beta-function ratio I,(p, q). Foster & Rees argue that the ‘generalized beta 
distribution’ is a (not the) natural generalization of the Beta distribution from uni- 
variate to multivariate analysis of variance; for other generalizations see [1], [2], 
[3], (4). 

The tables under review constitute a compilation of tables previously published 
in three papers by Foster and Rees [5], [6], [7]. Tabulated therein to 4D are values 
of the root of the equation 


I,(k; p,q) = P 
for P = .8(.05) .95, .99; 


k = 2, p = 3, 1(1) 10, g = 2 (1) (20) (5) 50, 60, 80; 
k = 3,4, p = 3(2) 4,q = 1 (1) 96. 


Two- to four-point Lagrangian interpolation in p and g is recommended; no specific 
accuracy is guaranteed. 

The computations for k = 2 were carried out on the N.R.D.C. Elliott 401 com- 
puter at Rothamsted; for k = 3, 4, on the DEUCE computer of the English Electric 
Company. Tables for P = .95, .99 and k = 2(1)6 have been given by Pillai {8}, [9]. 

Two examples are given of the application of the tables to the analysis of dis- 
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persion of means, and one example is given of such: application to the analysis of 
dispersion of regressions. 


J. ArtHuR GREENWOOD 


Princeton University 
Princeton, New Jersey 


. 8. S. Wixxs, ‘Certain generalizations in the analysis of variance,’’ Biometrika, v. 24, 
199, 'p. 471-494 
._ E. 8. Pearson & 8. 8. Wi1ks, ‘Methods of statistice' analysis appropriate for k samples 
of om variables,’’ Biometrika, v. 25, 1933, p. 353-378 


3. STATISTICAL RESEARCH Group, Cotumsra UNIVERsri x, Selected Techniques of Statistical 
— McGraw-Hill Book Co., New York, 1947, chap. 3, p. 111-184. 


. NerMaNn, Editor, Proceedings of the ‘Second Berkel y Symposium on M wiamationt Sta- 
titia and Probability, Univ. California Press, Berkeley & ios Angeles, 1951, p. 23 


5. F. G. Foster & D. H. Regs, ‘“‘Upper percentage poinis of the generalized beta. distribu- 
tion, 1,” Biometrika, v. 44, 1957, p. 237-247. [MT AC, Rev. -55, v. 12, 1958, p. 302] 
6. 


F. G. Foster, “Upper percentage points of the generalized beta distribution, IT,’ 
Biometrika, v. 44, 1957, p. 441-453 [[MTAC, Rev. 167, v. 12. 1958, p. 302] 


7. F. G. Foster, ‘ pper percentage — of the ge. atuined beta distribution, III,”’ 
ae v. 45, 1958, P. 492-503. [Math. Comp., Rev. 77, v. 14, 1960, p. 386] 
8. ‘0 


C.S8. Pruxar, mcise Tables for Statisticians, Statistical Center, Univ. of the Philip- 
pines, Mania 1957. 
9. 


K. C. 8. Pruuar & C. G. Banreaut, “On the distribution of the largest of six roots of a 
matrix in multivariate analysis,”’ Biometrika, v. 46, 1959, p. 237-240. 


43[K].—Tosto Krracawa & Micutwo Mrromg, Tables for the Design of Factorial 
Experiments, Dover Publications, Inc., New York, 1955 (printed in Japan; 
originally published by the Baifukan Company cf Japan as part 3 of the work 
with the same title), vii + 253 p., 26 em. Price $3.00. 


These tables consist of the actual tables that appeared in the original 1953 publi- 
cation in Japanese. An introduction to design principles and an explanation of the 
mathematical principles, parts 1 and 2 of the first p.ablication, have been omitted. 
Readers are now referred to Kitagawa’s Lectures on the Design of Experiments for 
this information and presumably for some help in the use of these tables. 

The American publisher’s jacket states that “this book contains tables for the 
design of factorial experiments and covers Latin squsres and cubes, factorial design, 
fractional replication in factorial design, factorial designs with split-plot confound- 
ing, factorial designs confounded in quasi-Latin squares, lattice designs, balanced 
incomplete block designs, and Youden’s squares.” Tie table of contents gives more 
detail under each of the eight main headings just listed, except for the last two. For 
example, orthogonal squares and cubes are listed, the 2” series of factorial arrange- 
ments goes up through 2°, mixtures of factorials such as a"b” mostly for m = 1 are 
listed, and the factorial replicates cover the 2” for 4, 3, }, and # replicates plus the 4 
replicate for 3". Perhaps it should be noted that tabies such as these are not really 
“for the design of experiments’’; the function of the tables is to help select a layout 
or make easy the randomization of the layout after the design has been selected. 

An examination of these tables shows that four Japanese pages have been cut 
out with scissors, and four English pages pasted in their place. The jacket further 
describes these tables as a ‘New revised edition. Explanatory notes.” The author’s 
preface does not describe what this reviewer would call a ‘revised edition’ and the 
explanatory notes consist of only one page. Since Kitagawa’s Lectures on Design 
of Experiments may not be readily available to some users of these tables, other 
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references might have been cited, e.g., O. Kempthorne, Design and Analysis of 
Experiments, W. G. Cochran and G. M. Cox, Experimental Designs, and O. L. 
Davies, Design and Analysis of Industrial Experiments. 

These tables are excellently and clearly printed. After one becomes acquainted 
with their structure and arrangement the tables should prove useful on many occa- 
sions to those persons engaged in the design of experiments in any field. One unique 
feature of these tables deserves notice. A complete listing of all 576 configurations 
of the 4x4 Latin square is given. Continuation of this procedure for larger squares 
would have produced a bulky volume. One wonders about the special utility of 
4x4 Latin squares which merited this complete listing. 

There are two comments that must be made about these tables. The first com- 
ment is a criticism on the failure to include a table of random numbers within the 
volume. This reviewer’s first act in using these tables will be to insert a small table 
of random numbers in both the front and rear of the volume. A table of experi- 
mental designs cannot be used without a random number table. As a consultant, 
when I pick up ‘my tables,’ I want to be sure that both items are with me. 

The second comment follows from the first. A preliminary section on randomiza- 
tion procedures and choice of specific layout for each design should have been in- 
cluded. If omitted, specific references to such instructions in the Fisher & Yates 
tables or in O. Kempthorne’s book should have been given. In this reviewer’s ex- 
perience both minor and major errors in designs have occurred because of a lack of 
clear understanding of proper randomization procedures. 

Finally, one may remark that these tables would have been much improved by 
the inclusion of some explanatory materials, and references for each design in- 
cluded. For statisticians, R. A. Fisher & F. Yates, Statistical Tables for Biological, 
Agricultural and Medical Research, Oliver & Boyd Ltd., Edinburgh (Fifth edition, 
1957), and E. S. Pearson & H. O. Hartley, Biometrika Tables for Statisticians, 
Vol. I, Cambridge, published for the Biometrika Trustees of the University Press 
(2nd printing, 1956) have set a high standard in this respect. The continued rapid 
development in the field of experimental design makes it difficult to keep tables of 
this type up to date. It is hoped that a really revised edition will soon appear. De- 
signs for response surface investigation and new fractional factorial arrangements 
need to be readily available. 


Emit H. Jese 


Operation Research Department 
Willow Run Laboratories 

The University of Michigan 
Ann Arbor, Michigan 


44[K].—Incram OLKIN, SupHisH G. Guurye, Wassitty Horrrpine, Wiii1AM G. 
Mapow, & Henry B. Mann, Editors, Contributions to Probability and Statistics, 
Essays in Honor of Harold Hotelling, Stanford University Press, 1960, x + 517 p., 
24 em. Price $6.50. ; 


This volume contains a collection of forty-two essays on probability and mathe- 
matical statistics in honor of Professor Harold Hotelling on his sixty-fifth birthday. 
The list of contributors, limited to those who have been closely associated with 
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Professor Hotelling, looks nevertheless like an up-to-date ““Who’s Who”’ in the sub- 
ject field. This fact alone pays an appropriate tribute to his influence and leadership. 

The first two essays, fitting to the occasion, deal with Hotelling the man, and 
as a leader and teacher in the field of mathematical statistics. The third one is a re- 
print of Hotelling’s own excellent paper on ‘“Teaching of Statistics,’ and the fourth 
one is a bibliography of his work. A total of ninety papers, not including reviews, 
were credited to him between 1925 and 1959—a truly impressive record of accom- 
plishment. 

The remaining thirty-eight research papers cover a wide spectrum of topics. 
There are seven papers on design and analysis of experiments, and about the same 
number in non-parametric statistics and also multivariate problems. Investigations 
into power, optimality, consistency, and robustness of tests, distribution theorems, 
and stochastic processes make up the bulk of the remaining papers. There is one 
paper on the inversion of partitioned matrices (Greenberg and Sarhan) and one on 
the numerical convergence of iterative processes (Moriguti). 

Since a listing of titles and authors takes about two pages, a detailed review of 
this diversified volume is an impossible task within the space allotted. If one paper 
has to be singled out as truly outstanding among the thirty-eight, I believe most 
people would agree to the choice of John Tukey’s ‘‘A Survey of Sampling from Con- 
taminated Distributions,” which investigates the robustness of efficiency of competi- 
tive estimators. In the paper the author considers two normal populations which 
have the same mean but whose standard deviations are in the ratio 3:1. One of the 
questions asked was: “‘What fraction of the wider normal population must be added 
to the narrower one in order for the mean deviation to be as good a large sample 
measure of scale as the standard deviation?” The answer, given two pages la:er, 
turns out to be a shockingly low .008. Tukey then suggests that ‘Problems of ro- 
bustness of efficiency are probably as important as problems of robustness of valid- 
ity, and, because of their relatively undeveloped stage, deserve even more attention 
from statisticians.’”’ No doubt this suggestion will be heeded. 

A list of titles and authors follows. Texts which are accompanied by tables are 
marked with an asterisk. The tables in paper No. 20 are separately described in the 
review immediately following. All the other tables are of illustrative nature, with 
limited selections of entries, and will not be discussed here. 


Part I. An Appreciation 
. Harold Hotelling—William G. Madow 
. Harold Hotelling—A Leader in Mathematical Statistics—Jerzy Neyman 
. The Teaching of Statistics—Harold Hotelling 
. Bibliography of Harold Hotelling 


Pe wn 


Part II: Contributions to Probability and Statistics 


5. Some Remarks on the Design and Analysis of Factorial Experiments—R. L. 
Anderson 

6. A Limitation of the Optimum Property of the Sequential Probabilit, Ratio 
Test—T. W. Anderson and Milton Friedman 

7. Decision Theory and the Choice of a Level of Significance for the t-Test— 
Kenneth J. Arrow 











ny 


nn 


ae a2 2h 2c Ae 








11. 
12. 


13. 
14. 


15. 
16. 


17. 
18. 


19. 


21. 
22. 


23. 
24. 
25. 
26. 
27. 


28. 
29. 
. Ranking in Triple Comparisons—R. N. Pendergrass and R. A. Bradley* 
31. 
32. 
33. 


34. 
35. 
36. 


37. 
. An Optimum Replicated Two-sample Test Using Ranks—Milton E. Terry* 
39. 
40. 
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. Simultaneous Comparison of the Optimum and Sign Tests of a Normal Mean— 


R. R. Bahadur 


. Some Stochastic Models in Ecology and Epidemiology—M. 8. Bartlett 
. Random Orderings and Stochastic Theories of Responses—H. D. Block and 


J. Marschak 

On a Method of Constructing Steiner’s Triple Systems—R. C. Bose* 

A Representation of Hotelling’s T* and Anderson’s Classification Statistic W in 
Terms of Simple Statistics—Albert H. Bowker 

Euler Squares—Kenneth A. Bush 

A Compromise Between Bias and Variance in the Use of Nonrepresentative 
Samples—Herman Chernoff 

Construction of Fractional Factorial Designs of the Mixed 2” 3” Series—W. 8. 
Connor 

Application of Boundary Theory to Sums of Independent Random Variables— 
J. L. Doob, J. L. Snell, and R. E. Williamson 

Some k-Sample Rank-order Tests—Meyer Dwass 

Characterization of Some Location and Scale Parameter Families of Distribu- 
tions—S. G. Ghurye 

Generalization of Some Results for Inversion of Partitioned Matrices—B. G. 
Greenberg and A. E. Sarhan 


. Selecting a Subset Containing the Best of Several Binomial Populations— 


Shanti S. Gupta and Milton Sobel* 

Consistency of Maximum Likelihood Estimation of Discrete Distributions—J. 
Hannan 

An Upper Bound for the Variance of Kendall’s ‘““Tau” and of Related Statis- 
tics—Wassily Hoeffding 

On the Amount of Information Contained in a o-Field—Gopinath Kallianpur 

The Evergreen Correlation Coefficient—M. G. Kendall 

Robust Tests for Equality of Variances—Howard Levene* 

Intrablock and Interblock Estimates—Henry B. Mann and M. V. Menon 

A Bivariate Chebyshev Inequality for Symmetric Convex Polygons—Albert 
W. Marshall and Ingram Olkin 

Notes on the Numerical Convergence of Iterative Processes—Sigeiti Moriguti 

Prediction in Future Samples—George E. Nicholson, Jr.* 


A Statistical Screening Problem—Herbert Robbins 

On the Power of Some Rank-order Two-sample Tests—Joan Raup Rosenblatt 

Some Non-parametric Analogs of ‘Normal’ ANOVA, MANOVA, and of 
Studies in ‘“Normal’”’ Association—S. N. Roy and V. P. Bhapkar 

Relations Between Certain Incomplete Block Designs—S. 8. Shrikhande 

Infinitesimal Renewal Processes— Walter L. Smith 

Classification Procedures Based on Dichotomous Response Vectors—Herbert 
Solomon ; 

Multiple Regression—Charles Stein 


A Survey of Sampling from Contaminated Distributions—John W. Tukey 
Multidimensional Statistical Scatter—S. S. Wilks 
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41. Convergence of the Empiric Distribution Function on Half-Spaces—J. Wolfo- 
witz 
42. Analysis of Two-factor Classifications With Respect to Life Tests—M. Zelen.* 
The five editors are to be congratulated for assembling and presenting this 
volume in an excellent manner. 
H. H. Ku 


National Bureau of Standards 
Washington, D. C. 


45[K].—Suant1 8. Gupta & Mitton SoskEt, “Selecting a subset containing the best 
of several binomial populations,” p. 224-248, Contributions to Probability and 
Statistics, Essays in Honor of Harold Hotelling, edited by Olkin et al., Stanford 
University Press, 1960. [See preceding review.] 


Given k binomial populations with unknown probabilities of success p: , p2 , 
-++ p, , a procedure R is studied by the authors which selects a subset that guaran- 
tees with preassigned probability P* that, regardless of the true unknown parameter 
values, it will include the best population; i.e., the one with the highest parameter 
value. Procedure FR for equal sample sizes is given as follows. Retain in the selected 
subset only those populations for which 7; = 2%max — d, where d = d(n, k, P*) is 
a non-negative integer, and x; denotes number of successes based on n observations 
from the ith population. Table 2 gives the values of d for k = 2(1)20, 20(5)50; 
n = 1(1)20, 20(5)50, 50(10)100, 100(25)200, 200(50)500; P* = .75, .90, .95, .99 
(a trial and error procedure R is given for large, unequal sample sizes) . 

Table 3 gives the expected proportion of populations retained in the selected sub- 
set by procedure R (for the special case pp = po = --: = Pia = DP, Mm = P+ 4, 
0<8<51,0S p <1 — 5) forn = 5(5)25; p* = .75, .90, .95; 5 = .00, .10, .25, 
50; and p + 6 = .50, .75, .95, 1.00. 

H. H. Ku 


National Bureau of Standards 
Washington, D. C. 


46[K].—Mavrice Henri QuENouILLE, “Tables of random observations from 
standard distributions,” Biometrika, v. 46, 1959, p. 178-202. Econ SHARPE 
Pearson, ‘‘Note on Mr. Quenouille’s Edgeworth Type A transformation,”’ 
Biometrika, v. 46, 1959, p. 203-204. 


Quenouille offers a random sample of 1000 each from the normal distribution 
and seven specified non-normal distributions. While a sample of 1000 is too small 
for much serious Monte Carlo work, the method of construction of the present 
tables, where the normal sample uniquely and monotonely determines the 7 non- 
normal samples, makes it suitable for pilot studies of the sensitivity of statistical 
procedures to departures from normality. 

Specifically, let 2, be a unit normal deviate from the tables of Wold [1]. Define 


ry 
y= (an) | exp (—42’) dz, 


ae = 3'7[2y — 1], 
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ts = 0.46271 e* — 0.76287, 
% —1 — log. (1 — ¥), 
124416 7, = —9552 + 127225 2, + 7824 2,” — 40 x, + 576 x, — 252 2’, 
M1 > —2.5, 


% = —1.86,% Ss —2.5, 
1536 a = 1411 2 + 562, — 32,’, 
124416 x; = —12144 + 122878 x, + 14304 x, — 1066 x,’ — 720 2,‘ + 261 z,’, 
a= —-2" log. [2 — 2y], u1>Od0 
as = 2" log, 2y, a = 0. 


Then, for i = 1(1)8, E(z;) = 0, E(x?) = 1. Here mis a rectangular random 
variate; x3, a log-normal variate; x, a one-tailed exponential variate; z,, a two- 
tailed exponential variate; 7; , x, 2; are Cornish-Fisher expansions with specified 
x; and «,. A short table on p. 179 shows that thespecifications are not met precisely; 
Pearson’s note shows that this failure is negligible for samples of 1000. 

The main table, p. 183-202, gives 1000 values of x;, i = 1(1)8, to 2 D, with 
Zz and =z’ in blocks of 50. Auxiliary tables on p. 180-182 give the first and second 
sample moments of the z; ; their theoretical x; , xs , xs , x¢ ; frequency distributions 
of the 8 samples; x; , 2%, x; to 3 D for x, = —3.2(.1) + 3.2. The italic headlines 
on p. 181-182 should be interchanged. 

It is not clear why random normal numbers were used as the basis for this table 
rather than random rectangular numbers, nor why the 2 D deviates of Wold [1] 
were chosen over the 3 D deviates of Rand Corp. [2]. 

Reprints may be purchased from the Biometrika Office, University College, 
London, W.C. 1, under the title “Tables of 1000 standardized random deviates from 
certain non-normal distributions.” Price: Two Shillings and Sixpence. Order New 
Statistical Tables, No. 27. 

J. ArtHuR GREENWOOD 


Princeton University 
Princeton, New Jersey 


1. Herman A. O. Wo.p, Random Normal Deviates. Tracts for Computers, no. 25, Cambridge 
Univ. Press, 1954. 
2. 


Tue RAND Corporation, A Million Random Digits With 100,000 Normal Deviates, The 
Free Press, Glencoe, Illinois, 1955. [MTAC, v. 10, 1956, p. 39-43). 


47(K].—ALrrep WeissperG & GLENN H. Beatty, Tables of Tolerance-Limit 
Factors for Normal Distributions, Battelle Memorial Institute, 1959, 42 p., 28 
cm. Available from the Battelle Publications Office, 505 King Avenue, Columbus 
1, Ohio. 


The abstract of the booklet reads as follows: “Tables of factors for use in com- 
puting two-sided tolerance limits for the normal distribution are presented. In con- 
trast to previous tabulations of the tolerance-limit factor K, we tabulate the factors 
r(N, P) and u(f, y), whose product is equal to K. This results in greatly increased 
compactness and flexibility. The mathematical development is discussed, including 
methods used to compute the tabulated values and a study of the accuracy of the 
basic approximation. A number of possible applications are discussed and examples 
given.” 

Since the mean yu and the standard deviation ¢ are frequently unknown, the toler- 
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ance limits must be computed on the basis of a sample estimate Z of the mean and an 
estimate s of the standard deviation. The tolerance limits treated in the booklet 
have the form x + Ks, where the factor K (the product of the tabulated entries 
r(N, P) and u(f, y)) accounts for sampling errors in ~ and s as well as for the 
population fraction P. 

Six levels of probability for P and y are used (.50, .75, .90, .95, .99, .999). The 
values of N used are given by 


N = 1(1)300(10)1000( 1000) 10000, ~. 
The values of f used are given by 
f = 1(1)1000(1000) 10000, ~. 


Values of r(N, P) and u(f, y) are given to four decimal places, which means that 
most of the tabular entries have five significant figures. 


Rosert E. GREENwWooD 


The University of Texas 
Austin, Texas 


48[K].—E. J. Wiuu1aMs, Regression Analysis, John Wiley & Sons, Inc., New York, 
1959, ix + 214 p., 24 em. Price $7.50. 


This useful volume is a monograph devoted to the exposition of the practical 
aspects of “regression analysis.”” These so-called regression analysis techniques are 
based on the method of least squares and are equivalent to analysis of variance pro- 
cedures. The author discusses many different techniques, some containing much 
novelty. All are accompanied by illustrations using actual data drawn mainly from 
the biological sciences. The book contains a great deal of interesting discussion and 
advice on the proper and practical applications of the methods. 

No attention is devoted to the planning of experiments; the book is only con- 
cerned with the analysis of data. Although nearly all the techniques involve the 
solution of simultaneous equations, there is little discussion of numerical tech- 
niques, except to recommend the “‘Crout’’ method. 

The author makes much use of statements about parameters which are termed 
fiducial statements. This reviewer feels these are confidence statements. In explain- 
ing the meaning of fiducial statements the author writes (p. 91), “...a fiducial 
statement about a parameter is, broadly speaking, a statement that the parameter 
lies in a certain range or takes a certain set of values. The statement is either true 
or false in any particular instance, but it is made according to a rule which ensures 
that such statements, when applied in repeated sampling, have a given probability 
(say 0.95 or 0.99) of being correct.” 

The various techniques are presented without theory, as ‘“‘to have done so would 
have made the book unnecessarily long.’’ Without the accompanying theoretical 
material, this book is simply a handbook of regression methods. It is for this reason 
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that the reviewer feels the book will be more useful to applied statisticians than to 
the author’s intended audience, i.e., research workers in the experimental sciences. 


M. ZELEN 


University of Maryland 
College Park, Maryland 


49[L].—G. F. Miter, Tables of Generalized Exponential Integrals, National Physi- 
cal Laboratory Mathematical Tables, Vol. 3, British Information Services, New 
York, 1960, iii + 43 p., 28 cm. Price $1.43 postpaid. 


According to the author, the tables under review were prepared to meet the 
requirements of quantum chemists concerned with the evaluation of molecular 
integrals, who frequently have found the tables computed by the New York Mathe- 
matical Tables Project and edited by G. Placzek [1] inadequate for this purpose. 

Actually tabulated in the present work is the auxiliary function 


F(z) = (x + n)eE,(z), 
where E,(z) represents the generalized exponential integral, defined by the equa- 
tion E,(2z) -/ e 't" dt. Three tables of F,,(z) to 8D are provided. The first table 
1 


covers the range x = 0(0.01)1 for nm = 1(1)8; the second, the range z = 0(0.1)20 
for n = 1(1)24; and the last, the range 1/x = 0(0.001)0.05 for n = 1(1)24. Modi- 
fied second (and occasionally fourth) central differences are provided throughout 
for use with Everett’s interpolation formula. For details of methods of interpolation 
and tables of interpolation coefficients the table-user is referred to the first two 
volumes of this series of tables [2], [3]. 

It is stated that the total error in an unrounded interpolated value of F,(x) 
derived from the present tables need never exceed 14 units in the eighth decimal 
place if the tabulated differences are used. Furthermore, the values of F(z) here 
tabulated are guaranteed to be accurate to within 0.6 unit in the last place. 

The tables are preceded by an Introduction containing a brief account of perti- 
nent literature, followed by a section devoted to a description of the tables and a 
justification for the tabulation of F(z) in preference to E,(z). The properties of 
the generalized exponential integral, many of them reproduced from Placzek [1], 
are enumerated in a third section. The fourth section of the text is devoted to a 
careful description of the several procedures followed in the preparation of the 
tables. An excellent set of references is appended to this introductory textual 
material. 

The typography is uniformly excellent, and the format of the tables is conducive 
to their easy use. The only defect observed was a systematic error in the heading 
of Table 2 on pages 24 through 37, where this heading erroneously appears as 
Table 3. 

J. W. W. 


1. NatronaL Researcu Councit or Canapa, Division of Atomic Energy, Report MT-1, 


The Functions E,(x) -/ e~=“y-" du , Chalk River, Ontario, December 1946. Reproduced in 
1 
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Nat. Bur. Standards A pl Math. Ser. No. 37, Tables of Functions and of Zeros of Functions, 
— a See RM 392, MT AC, v. 2, 1946-47, p. 272; and RMT 104, MT AC, v. 10, 1956, 
p. 249— 


2. L. Fox, The Use and Construction of Mathematical Tables, National Physical Laboratory 
Mathematical Tables, v. 1, London, 1956. See RMT 8, MTAC, v. 13, ae: 61-64. 
Aa 


3. L. Fox, Tables of Everett Interpolation Coefficients, National Physical Laboratory Mathe- 
matical Tables, v. 2, London, 1958. 


50[L].—-F. W. J. Otver, Editor, Bessel Functions, Part III, Zeros and Associated 
Values, Royal Society Mathematical Tables No. 7, Cambridge University Press, 
New York, 1960, lx + 79 p., 29 em. Price $9.50. 


The present volume is a step towards the completion of a program for the tabu- 
lation of Bessel functions initiated by the British Association Mathematical “Tables 
Committee, and continued since 1948 by the Royal Society Mathematical Tables 
Committee. Part I of this series, Bessel Functions, Functions of Order Zero and Unity 
appeared in 1937, and Part II, Bessel Functions, Functions of Positive Integer Order 
appeared in 1952 (see MTAC v. 7, 1953, p. 97-98). Recall that Part I contains a 
section on the zeros of J,(z), Yn(z), m = 0, 1, but Part II is without a section 
devoted to zeros. 

Part III, the present work, deals with the evaluation of zeros of the Bessel func- 
tions J,(z) and Y,(z) for general » and z. Tables are also provided as described 
later in this review. A history of the project is given in the “Introduction and 
Acknowledgements,”’ by C. W. Jones and F. W. J. Olver. A chapter on “Definitions, 
Formulae and Methods” by the above authors is a valuable compendium of tech- 
niques for the enumeration of zeros and associated functions. In particular, it is an 
excellent guide if zeros are required of other transcendental functions which satisfy 
second-order linear differential equations. Several methods of computation are out- 
lined. For instance, the method of McMahon is useful for v fixed and z large, while 
the inverse interpolation approach of Miller and Jones presupposes a tabulation of 
the functions themselves. Between the regions covered by these techniques is a gap 
which increases with increasing v. The gap is bridged by application of Olver’s 
important contributions on uniform asymptotic expansions of Bessel functions. 

The section ‘Description of the Tables, Their Use and Preparation”’ is by the 
editor. A short description of the tables follows. Table I gives zeros j,,, of Jn(x), 
Yn, of Y,(x), and the values of Jn’(jn,2), Yn'(Yn.s). Table II gives zeros j,,, of 
J,'(x), Yn. Of Y»'(x), and the values of Jn(jn.2), Yn(Yn.s). Table III gives zeros 
Am, , Um,» Of the derivatives jm’(x), ym’(x) of the spherical Bessel functions jn(x) = 
(x/2a)"? SF mare(t), Ym(t) = (4/2x)"?¥ngie(x), and the values of jim(Om,s) 
Ym( bm»). The ranges covered are 


0(4)203,  s = 1(1)50, Tables I and II; 
0(1)20,  s=1(1)50, Table III. 


n 


m 


All entries are to eight decimals, and in no case should the end-figure error exceed 
0.55 of a unit in the eighth decimal. 

The coefficients in the uniform asymptotic expansions (previously mentioned ) 
which are used to evaluate items in Tables I-III for n (or m) large are given in 
Table IV. The expansions for the Bessel functions of Tables I-III also depend on 
zeros and associated values of certain Airy functions and their derivatives. These 
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data are recorded in Table V. Further description of the contents of ‘Lables IV-V 
requires much more space and so is omitted here. Suffice it to say that with the 
aid of these tables, the entries in Tables I-III can with a few exceptions be evalu- 
ated to at least eight significant figures for 20 <n << ~,lSs<@ 

There is a good set of references. The printing and typography are excellent, 
and the present volume upholds the eminent tradition of British table-makers. 


Y. L. L. 


51[L, V).—J. W. Mixes, “The hydrodynamic stability of a thin film of liquid in 
uniform shearing motion,” J. Fluid Mech. 8, Pt. 4, 1960, p. 593-610. (Tables 
were computed by David Giedt.) 

Let 


§(z) = [1 — F(z)|" = w[As(—w)]" E T [ A,(-t) a|, w = ze™"®, 


F(z) = 2 F(z) + we [A,’(—w)] "A —w) [1 — F(z)]. 
6 (z) = (2) + F(z), k=0,1; F(z) = Pr(z) + iF((z). 


The paper contains tables of §(z), ¥’(z), F(z) and 2°F;(z) for z = —6(.1)10, 48. 
The tables were obtained on an automatic computer by numerical integration of 
an appropriate differential equation. It can be seen from the above that the tables 
depend on values of the Airy integral A,(z), its derivative and integral along the 
rays +/6 and —5z/6 in the complex plane. Tables of A,(z) and its derivative are 
now available for complex z in rectangular form, but not in polar form. Also, tables 


of [ A,(-+t) dt are available for z real. Thus, the given tables depend on values of 
0 


some basic functions which, if available, would cut new ground. Unfortunately, the 
basic items were swallowed up in the automatic computation of F(z). We have here 
a poor example of table making,—a practice which should not be emulated. 


Y. L. L. 


52[P].—Hetmut Horses (Compiler), Wasserdampftafel der Allgemeinen Elektricitats- 
Gesellschaft, R. Oldenbourg, Munich, 1960, 48 p., 30 em. DM 16 (Paperback). 


There are two tables in this collection. Table I is a four-place table giving the 
temperature, the specific volume, the specific enthalpy, and the specific entropy as 
functions of the absolute pressure p. The last three dependent variables are given 
both for the fluid state and the gaseous state. The variable p ranges from 0.010 to 
225,650 atmospheres, and the interval varies from 0.001 to 2000. Table II gives 
the specific volume, specific enthalpy and specific entropy as functions of tempera- 
ture for constant pressure. Here p has the values 1, 5, 10 (10) to 400 atmospheres, 
and ¢ varies from 0 (10) to 330 degrees centigrade. 

The tables were calculated by expressing each of the dependent variables as 
polynomials in the pressure with coefficients as functions of the temperature or in 
some cases functions of the temperature and pressure. The error bounds given by 
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such approximations to the specific volume and specific enthalpy as functions of 
pressure and temperature are included in the collection. 


A. H. T. 


53[P, X, Z].—Warp C. Saneren, Digital Computers and Nuclear Reactor Calcula- 
tions, John Wiley & Sons, New York, 1960, xi + 208 p., 24 em. Price $8.50. 


As the author states in his preface, the primary objective of this book is to pre- 
sent to nuclear engineers and scientists an introduction to high speed reactor calcu- 
lations. Since the appearance of the basic reference, The Elements of Nuclear Reactor 
Theory by Glasstone and Edlund, Van Nostrand, 1952, the entire complexion of 
actual reactor design calculations has changed as a result of the growth in speed and 
size of computing machines, and reactor design calculations represent today a sig- 
nificant part of scientific computing time on modern computers. 

The outline of the book by chapters is 

Chapter 1. Introduction 

Chapter 2. Digital Computers 

Chapter 3. Programming 

Chapter 4. Numerical Analysis 

Chapter 5. A Code for Fission-Product Poisoning 
Chapter 6. Diffusion and Age-Diffusion Calculations 
Chapter 7. Transport Equation—Monte Carlo 
Chapter 8. Additional Reactor Caiculations 

In Chapter 1, the author reviews the tremendous parallel growth of high speed 
computing machines and nuclear reactors, and their present interplay. In Chapter 
2, an introduction and description of present day computers is given. In Chapter 3, 
programming for computers is introduced. After some preliminary remarks (no 
proofs) about interpolation, numerical integration, matrices, etc., items which can 
be found in many well-known texts on elementary numerical analysis, the author 
treats in Chapter 4 the more relevant problem of the numerical approximation of 
partial differential equations by difference equations, and their solution by means 
of iterative methods. Also, the treatment of interface conditions, which arise 
naturally in heterogeneous reactors, is given. 

In Chapter 5, a simple code for fission-product poisoning is followed from the 
physical and mathematical definitions through to the construction of a program in 
the Bell (Wolontis) system. 

In Chapter 6, the longest chapter, the author describes diffusion calculations, 
extending from steady-state criticality problems for reactors to the solution of two- 
and three-dimensional multigroup diffusion equations. In Chapter 7, the S,, method 
of Carlson is described, along with the use of Monte Carlo methods for solving prob- 
lems such as those encountered in shielding calculations. 

In his primary aim, the author does succeed. Nevertheless, the reviewer, being 
quite familiar with this area, was most critical with respect to the age of the refer- 
ences, as most of the technical papers referred to had appeared prior to 1957. As 
no serious attempt was made to fill the gap between these earlier developments and 
the developments which have taken place in the reactor field in the last few years, 
many statements in the book are either somewhat obsolete or misleading. For 
example, the numerical inversion of tridiagonal matrix equations on page 74 by an 
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algorithm is not stated to be simply Gauss elimination applied to the matrix prob- 
lem, and in fact the author states that this “method” has not appeared in textbooks 
as yet. The iterative methods of Young-Frankel, and Peaceman-Rachford are each 
discussed twice, (p. 84 and p. 144) and not one of the four definitions is completely 
accurate. The book is, however, the only existing bridge between The Elements of 
Nuclear Reactor Theory and present computational technique in the reactor field. 


R. 8. V. 


54[S, W].—Ronatp A. Howarp, Dynamic Programming and Markov Processes, 
Technology Press & Wiley, New York, viii + 136 p., 23 em. Price $5.75. 


Consider a physical system S represented at any time t by a state vector z(t). 
The classical description of the unfolding of the system over time uses an equation 
of the form z(t) = F(2(s), s = t), where F is a prescribed operation upon the 
function z(s) for s S ¢t. In certain simple cases, this reduces to the usual vector 
differential equation dz/dt = g(x), x(0) = c. 

For a variety of reasons, it is sometimes preferable to renounce a deterministic 
description and to introduce stochastic variables. If we take z(t) to be a vector 
whose i-th component is now the probability that the system is in state i at time t, 
and allow only discrete values of time, we can in many cases describe the behavior 
of the system over time quite simply by means of the equation z(t + 1) = Az(t). 
Here A = (a;;),7,j7 = 1,2, --- , N, is a transition matrix whose element a,; is the 
probability that a system in state j at time ¢ will be found in state i at time ¢ + 1. 
Processes of this type are called Markov processes and are fundamental in modern 
mathematical physics. 

So far we have assumed that the observer plays no role in the process. Let us 
now assume that in some fashion or other the observer has the power to choose the 
transition matrix A at each stage of the process. We call a process of this type a 
Markovian decision process. It is a special, and quite important, type of dynamic 
programming process; cf. Chapter XI of R. Bellman, Dynamic Programming, 
Princeton University Press, 1957. 

Let us suppose that at any stage of the process, we have a choice of one of a 
set of matrices, A(q) = (a;;(q)). Associated with each choice of g and initial state i 
is an expected single-stage return b;(q). We wish to determine a sequence of choices 
which will maximize the expected return from n stages of the process. Denoting 
the maximum expected return from an n-stage process by f;(n), the principle of 
optimality yields the functional equation 


fi(n) = rn [bi(q) + 2; assaf — 1). 


In this form, the determination of optimal policies and the maximum returns 
is easily accomplished by means of digital computers; see, for example S. Dreyfus, 
J. Oper. Soc. of Great Britain, 1958. Problems leading to similar equations, resolved 
in similar fashion, arise in the study of equipment replacement and in continuous 
form in the “optimal inventory” problem; see Chapter Five of the book mentioned 
above and K. D. Arrow, 8. Karlin, and H. Scarf, Siudies in the Mathematical Theory 
of Inventory and Production, St.nford University Press, 1959. 
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As in the case of the ordinary Markov process, a question of great significance 
is that of determining the asymptotic behavior as n — ~. It is reasonable to suspect, 
from the nature of the underlying decision process, that a certain steady-state 
behavior exists as n — «©. This can be established in a number of cases. 

The author does not discuss these matters at all. This is unfortunate, since there 
is little value to steady-state analysis unless one shows that the dynamic process 
asymptotically approaches the steady-state process as the length of the processes 
increases. Furthermore, it is essential to indicate the rate of approach. 

The author sets himself the task of determining steady-state policies under the 
assumption of their existence. Granted the existence of a “steady state,”’ the func- 
tions f;(n) have the asymptotic form nc + b; + 0(1) asn — o, where c is inde- 
pendent of 7. The recurrence relations then yield a system of equations for c and 
the b; - 

This system can be studied by means of linear programming as a number of 
authors have realized; see, for example, A. S. Manne, ‘“‘Linear programming and 
sequential decisions,’ Management Science, vol. 6, 1960, p. 259-268. 

Howard uses a different technique based upon the method of successive approxi- 
mations, in this case an approximation in policy space. It is a very effective tech- 
nique, as the author shows, by means of a number of interesting examples drawn 
from questions of the routing of taxicabs, the auto replacement problem, and the 
managing of a baseball team. 

The book is well-written and attractively printed. It is heartily recommended 
for anyone interested in the fields of operations research, mathematical economics, 
or in the mathematical theory of Markov processes. 

RicHARD BELLMAN 


The RAND Corporation 
Santa Monica, California 
(Reprinted from Quarterly of Applied Mathematics) 


55[V, X].—Gerrrupe Biancn, Kart Gorrrriep GuDERLEY & Emma MARIAN 
VALENTINE, Tables Related to Axial Symmetric Transonic Flow Patterns, WADC 
Technical Report 59-710, 1960, Office of Technical Services, U. 8. Department 
of Commerce, Washington 25, D. C., xlviii + 108 p., 27 cm. 


The equations of motion of a compressible fluid are non-linear and are generally 
difficult to handle. In certain cases, such as in the flow past slender bodies of revo- 
lution, the equations can be approximated by much simpler ones. For subsonic and 
supersonic flow these approximating equations are linear. When the flow velocity 
is nearly equal to one, the approximate equation for the disturbance potential takes 
the non-linear form 


—, &,, + %,, 4% + 3 D..5 = 0 
y' ¥ 


when x, y, w are cylindrical coordinates. K. G. Guderley and his colleague H. 
Yoshihara have studied the flow past slender bodies at Mach numbers close to one 
in a series of papers and in a book by Guderley, Theorie schallnaher Strémungen, 
Springer-Verlag, 1957. 
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The basic technique applied to axially symmetric flows is to find a basic solution 
©” of the form 


@” = y'" F(r) 


where ¢ = x/y” and n is a constant. The variable f then satisfies the ordinary dif- 
ferential equation 


(f’ — n't?)f” + (5n? — 4n)tf’ — (3n — 2)f = 0. 


Further solutions necessary to satisfy particular boundary conditions are then 
found by perturbing the basic solution by the function 4, i.e., 


& = + Bz, y, w). 


Then @ is assumed to satisfy the linear equation 


“i - - & . b, 
—,"6,, = : o, + %,, + = + y¥ = 0. 


Particular solutions of the equation are then found in the form 
& = y'g(f) cos ma; m = 0,1, 2,--- 
Then g({) satisfies the ordinary differential equation 
(f’ — n'g*)g” + [f” + (2un — n*)e\g’ + (m’ — v*)g = 0. 


The solution of this equation leads to an eigenvalue problem with the eigenvalue v. 

The present report tabulates the functions f and g together with their deriva- 
tives and some other related functions. In Table 1 appear 6D values of f(¢) and df/dt 
for ¢ = —7.5(.1) — 3(.02)1; in Table 2 similar information appears for g(¢) and 
dg/d¢. These tabular data are given for sever} values of the eigenvalue ». 

A rather complete discussion of the mathematical problems involved is given 
in the introduction. The eigenvalues are found using a contour integration tech- 
nique. It is stated that the numerical calculations are performed on an ERA 1103, 
with considerable pains taken to insure accuracy. The entries are stated to be 
correct to within one unit in the last place. 

The tables should be quite useful to anyone interested in the study of special 
cases of transonic flow. 

Ricuarp C. RoBerts 


U. 8. Naval Ordnance Laboratory 
White Oak, Silver Spring, Maryland 


56[W].—H. L. Tooruman, A Table of Probability Distributions useful in War 
Games and other Competitive Situations, NRL Report 5480, U.S. Naval Research 
Laboratory, Washington, D. C., May 16, 1960, i + 91 p., 27 em. 


A player makes a maximum of (2r — 1) plays. On odd-numbered plays he 
scores | with probability p, and 0 with probability 1 — p, ; on even-numbered plays 
he is eliminated from subsequent play with probability p:.. The probability that 
he will score exactly n is 


r—1 
nf rn r— k 3 —s — 
S, = (") pr"(1 — pr)” "(1 — po)” + > (F) ata — pr)” "pall — pr)*. 
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The probability, U,, that m independent players score a total of exactly ¢ is the 
coefficient of x‘ in (2S,2")”. 

The table on p. 6-91 gives U, to 4D for m = 1(1)4,r = 1(1)4,¢ = O(1)mr. 
For p; = .01(.01).06(.02).22, .25, p, = 0(.05).2(.1).9; for p: = .3(.05).95, p. = 
0(.01).02(.02).12, .15(.05).9. U, was computed to 9D or better on the NAREC, 
and each value was rounded to 4D individually; i.e., the U; were not forced to sum 
to 1. Quadratic interpolation in p; or pe is stated to yield a maximum error of .0016. 

The typography (photo-offset reproduction of Flexowriter script) is adequate 
but undistinguished; all decimal points are omitted from the body of the table. 


J. ARTHUR GREENWOOD 


Princeton University 
Princeton, New Jersey 


57[W, X].—S. Vaspa, An Introduction to Linear Programming and the Theory of 
Games, John Wiley & Sons, Inc., New York, 1960, 76 p., 22 em. Price $2.25. 
This book introduces the basic mathematical ideas of linear programming and 

game theory (mostly matrix games) in a form suitable for anyone who has had a 
little analytic geometry (and is not frightened by subscripts and double subscripts). 
Part I, on linear programming, begins with two examples, the second of which is a 
transportation problem, and then describes the simplex method of solving the 
transportation problem. Then comes the graphical representation of the general 
linear programming problem, followed by the general simplex method and a dis- 
cussion of such complications as finding a first feasible solution, multiple solutions, 
and degeneracy. The chapter closes with the duality theorem. 

Part II, on games, begins with two examples of matrix games, the second of 
which admits no saddle point, and introduces the concepts of mixed strategy and 
value. This is followed by a discussion of games in extensive form, and their normali- 
zation. A section on graphical representation is followed by the description of the 
equivalent linear program, and the Shapley-Snow “algorithm” is offered as an 
alternative method of calculating equilibrium strategies and value. Next the concept 
of equilibrium point in non-zero sum games is discussed, followed by three examples 
of infinite games. The book closes with an appendix proving the main theorem of 
matrix games along Ville’s lines. 

This book, compiled from lecture notes of short courses offered by the author, is 
suitable as a text for a short course for students with slight mathematical prepara- 
tion. 


ALAN J. HorrmMan 


General Electric Co. 
New York 22, N. Y. 


58[W, Z].—Franz L. Aur, Editor, Advances in Computers, Vol. 1, Academic Press, 
Ine., New York, 1960, x + 316 p., 24 em. Price $10.00. 


Advances in Computers is a useful addition to the rapidly growing literature on 
modern high-speed computers and their application. It is intended by the editor to 
be the first volume in a series which will contain monographs by specialists in vari- 
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ous areas of work in this field. These are to be written in non-technical language, so 
as to be easily understood by specialists in areas other than those of the writer. 

The present volume contains six articles by well-qualified authors in six signifi- 
cant and interesting areas of work related to computers. Four summarize progress 
to date in the application of computers to weather prediction, translation of lan- 
guages, playing games, and recognition of spoken words. Two are related to tech- 
niques used in computer programming and design. The titles and authors are: 

1. General-Purpose Programming for Business Application—Calvin E. Gotlieb 

2. Numerical Weather Prediction—Norman A. Phillips 

3. The Present Status of Automatic Translation of Languages—Yehoshua Bar- 
Hillel 

4. Programming Computers to Play Games—Arthur L. Samuel 

5. Machine Recognition of Spoken Words—Richard Fatehchand 

6. Binary Arithmetic—George W. Reitwiesner. 

Since most of the areas of work covered by the papers in this volume are in a 
rapid state of flux, the assignment to write survey papers in these areas, undertaken 
by the authors, is a most difficult one. Each author has proceeded to carry out this 
assignment in his own characteristic manner. Thus, Gotlieb attempts to present a 
factual summary of some of the programming procedures used at present in process- 
ing data for business applications; whereas, Yehoshua Bar-Hillel presents a critical 
evaluation of the various efforts conducted in the field of automatic translation of 
languages—at times, highly critical. A large part of the material covered is ad- 
mittedly subjective, and bears the imprint of the writers’ points of view and contri- 
butions. Nevertheless, the six papers in this volume constitute authoritative sur- 
veys of the areas of work discussed. Together with the bibliographies given at the 
end of each paper, these articles will be valuable to the new researcher in the fields 
covered, as well as to the interested layman who wishes to familiarize himself with 
the exciting advances in computer technology. 

H. P. 


59[Z].—AnpreEw D. Bootu, Automation and Computing, The Macmillan Co., New 
York, 1959, 158 p., 21 em. Price $5.00. 


This book is intended mainly for the educated layman. In it the author attempts 
“to bridge the gap between the superficial accounts of electronic computers and 
automation ...and the specialists’ monographs. ...”’ He has given an admirably 
written and lucid account of digital and analogue computers. His three chapters on 
the logical design of digital computers, the physical basis of this design, and pro- 
gramming for digital computers are very clear and informative, though concise. 

The three chapters on automation in clerical work, control of continuous proc- 
esses, and automatic machine tools and assembly processes are not as well done as 
the first three. The well-educated layman will have to expend a great deal of effort 
in order to follow the discussion in these ehapters. 

The last two chapters entitled ‘Strategic and Economic Planning” and “‘Non- 
numerical Applications of Computing Machines” are very brief. The former is much 
too short to give the reader more than a glimmer of what is involved in game 
theory. The last chapter furnishes a well-written introduction to methods for non- 
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numerical applications of computers, but, because it is so short, leaves the reader 
wishing the author had devoted more space to this subject. This reader would have 
preferred to have the author do this and omit some of his pronouncements on gov- 
ernment (for example, the discussion on page 20 beginning with “Democratic govern- 
ment, too, is an example of Man in decay, .. .’’). 

There are a few typographical errors in the book. The most disturbing one ap- 
pears on page 36 where the binary addition table has the entry 


1 + 1 = 1 (carry 1). 
AT. 


60[Z].—Rosert H. Gregory & Ricnarp L. Van Horn, Automatic Data-Processing 
Systems, Wadsworth Publishing Co., San Francisco, 1960, xii + 705 p., 23 em. 
Price $11.65. 


This introductory book on automatic data-processing systems (ADPS) is a re- 
vision of a text which was used in management development courses sponsored by 
the Army Ordnance Corps. The affirmative objective is to instruct, enlighten, and 
inform management on the developments, techniques and applications of methods 
in management science, mathematics, and large-scale computing for the solution of 
today’s complex business problems. 

The book is divided into seven parts and three appendices. In Part One, ‘‘Orien- 
tation,” the principles of basic computer programming are elucidated by means of 
a hypothetical computer which embodies an instruction repertoire of several exist- 
ing machines. Various numerical and alphanumerical coding systems for storing 
data on punched cards, punched paper tapes, and magnetic tapes are also discussed 
here. 

Part Two, “Automatic Equipment,” deals with input-output hardware, storage 
devices, arithmetic and control units. The section concludes with a synopsis of the 
salient characteristics of approximately twenty computing systems: speed, storage, 
instruction repertoires, tapes, and peripheral equipment. 

Advanced programming techniques and systems provide the subject matter of 
Part Three, “Programming and Processing Procedures.” In this section the authors 
present a synthesis of the pros and cons of automatic programming and integrated 
data processing, two important and topical subjects. 

The role of the data-processing unit in ‘management information systems’’ is 
the theme of Part Four, ‘Principles of Processing Systems.’”’ Several methods are 
suggested for selecting from a welter of available facts the pertinent information for 
effective executive decision-making. The reporting-by-exception principle is de- 
scribed in detail. Since the efficacy of the final system design is inextricably related 
to economic considerations, the authors analyze the major factors for determining 
the cost of obtaining and processing data, and explore the concept of the ‘“‘value”’ of 
information in relation to its cost. The last chapter in Part Four outlines the broad 
principles underlying systems analysis and design. 

Factors that affect the organizational structure of data processing are subject to 
examination in Part Five, “Systems Design.” In particular, considerable attention 
is devoted to problems associated with centralized data processing and decentralized 
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management control. General tools for systems analysis and specific data-processing 
techniques are also described here. 

Part Six, ‘Equipment Acquisition and Utilization,” presents in a nontechnical 
manner a methodical approach for evaluating, selecting, installing, and implement- 
ing automatic data-processing systems for business problems. Considerable space is 
devoted to the preparation of feasibility studies, application studies, and equipment 
acquisition proposals. This is followed by a detailed discussion of the problems en- 
tailed in the installation of new equipment. 

The concluding portion, Part Seven, “System Re-examination and Prospective 
Developments,” touches on a variety of mathematical techniques for the solution 
of management problems and concludes with a discussion on anticipated future 
developments in automatic data processing. 

Three appendices are: 

I. History of Computation and Data-Processing Devices 
II. Questions and Problems 

III. Glossary of Automatic Data-Processing Terminology 

Although the treatment of the basic principles of computer programming illu- 
minates many of the complex and important aspects of business data processing, 
the authors give little heed to the practical requirements for large-scale production 
system runs. Such concepts and techniques as rerun procedures, interior tape labels, 
alternation of servos, and programming methods for effective utilization of buffer- 
ing are not even mentioned, while the subjects of editing, flow charting of instruc- 
tion routines, and sorting techniques for large tape files are glossed over. But on the 
whole, the informative and lucid presentation of the general principles of automatic 
data processing from the standpoint of business systems will provide management 
personnel with a short, intensive, and enlightened education on electronic com- 
puters and their impact on business data processing. 

MILTON SIEGEL 
Applied Mathematics Laboratory 


David Taylor Model Basin 
Washington 7, D. C. 


61[Z].—Antuony G. Ortrincer, Automatic Language Translation, Harvard Uni- 
versity Press, Cambridge, 1960, xix + 380 p., 24 cm. Price $10.00. 


Automatic Language Translation by Anthony G. Oettinger is the eighth in a 
series of Harvard Monographs in Applied Science. ‘“These monographs are devoted 
primarily to reports of University research in the applied physical sciences, with 
special emphasis on topics that involve intellectual borrowing among the academic 
disciplines.” Professor Oettinger’s monograph is devoted to the lexical and technical 
aspects of automatic language translation, with particular emphasis on Russian-to- 
English translation. 

The contents of this work can quickly be conveyed by the titles of its chapters. 
Chapter 1, “Automatic Information-Processing Machines,” discusses the organiza- 
tion, elements of programming, and the characteristics of information-storage 
media. Chapter 2, ‘The Structure of Signs,” differentiates the notions of use, men- 
tion, and representation of signs; mathematical transformations; and mathematical 
models. Chapter 3, ‘““Flow Charts and Automatic Coding,” treats the use of flow 
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charts, addressing, algorithms, and programs. Chapter 4, ““The Problem of Trans- 
lation,” is a general discussion of the possibilities and types of automatic transla- 
tion, grammar, and interlingual correspondence, and syntactic and semantic prob- 
lems. Chapter 5, “Entry Keys for an Automatic Dictionary,” discusses inflection, 
paradigms, affixes, stems, and inflection algorithms for Russian. Chapter 6, ‘‘Mor- 
phological and Functional Classification of Russian Words,” consists of a detailed 
account of Oettinger’s morphological and functional system for Russian words; 
nominal forms, adjectival forms, and verbal forms; and an appendix that gives synop- 
tic classification tables. Chapter 7, “Dictionary Compilation,” describes the 
structure of the Harvard automatic dictionary by giving the structure of items and 
files, methods of detecting and correcting mistakes in transcription and classification, 
and English correspondence and grammatical codes. Chapter 8, ‘Dictionary Opera- 
tion,” describes the function of the Harvard automatic dictionary, lookup pro- 
cedures, and word-by-word translation ; and is followed by an appendix that presents 
an edited trot, the transcription of the edited trot, and an example of conventional 
translation. Chapter 9, “Problems in Dictionary Compilation and Operation,” dis- 
cusses the problems of paradigm homography, stem homography, “short’’ words, 
and a detection and correction of mistakes in dictionary compilation. Chapter 10, 
“From Automatic Dictionary to Automatic Translator,” presents the author’s 
views on how the Harvard automatic dictionary might lead to a complete system of 
automatic translation. 

Since this is the first book published in America devoted to automatic translation 
of languages, ‘i is a landmark. Several cautions should be mentioned, however, for 
those who are not familiar with the state of progress in machine translation. First, 
this book is not a work devoted to the general problem of translating one natural 
language to another. It is highly specialized, since it treats only the Russian-to- 
English translation problem. Second, much of the book is devoted to the very 
detailed description of the particular automatic dictionary compiled at Harvard 
University. This description does not permit conclusions to be drawn “auto- 
matically” about dictionary compilation at other machine translation research 
centers. Third, all detailed computer descriptions are in terms of the Sperry-Rand 
UNIVAC I computer, whereas almost all other machine translation programs in 
the United States are written for IBM 704 or 709 computers. 

Nevertheless, Professor Oettinger is to be congratulated for presenting the first 
detailed, and scientifically accurate description of any machine translation project 
in the U.S., if not in the world. As such, this book will be of interest. to computer 
scientists, mathematicians, linguists, and to others interested in acquiring knowl- 
edge about this important subfield of modern linguistic analysis. As this reviewer 
has often emphasized, the gains to be made in linguistic analysis will overshadow 
those which have been made in numerical analysis. 

H. P. EpmMunpDsoNn 


Planning Research Corporation 
Los Angeles 24, California 
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TABLE ERRATA 


299.—Mme. JacqueLine Heurtaux, “Tables de polynémes d’interpolation avec 
seulement deux abscisses distinctes,” Chiffres, 1" Année, Paris, March 1958, p. 
25-34. 


) for read 
p. 31, Qe, 2=0.15 0.97338 82250 0.97338 81250 
H. E. Sauzer 


Convair-Astronautics 
San Diego 12, California 


300.—H. Taxryama, “Expressions for interpolation and numerical integration of 

high accuracy,” Tohoku Univ. Technol. Reports, vy. XXIII, 1958, p. 47-70. 

On p. 69, corresponding to u = 0.04, the value of U,’ should read 0.039... 
instead of 9.039... ; and corresponding to u = 0.34, the value of U, should read 
0.72203 53338 6336 instead of 0.72203 58338 6336. 

H. E. Sauzer 


CORRIGENDA 


C. W. Dunnett & R. A. Lamm, “Some tables of the multivariate normal probability 
integral with correlation coefficients 4,’’ Math. Comp., Review 50, v. 14, 1960, 
p. 290. 


In the expression given for the probability integral of the multivariate normal 
distribution in n dimensions the upper limit of the innermost integral should read 


x, instead of x», and the denominator (1 — p)“>~ should be replaced by 
{a=1) 
(i—~») ° 
In the following line of the text 
jor Fugli,°**, a), eal Fiala, *** . Pu). 


F. R. GANTMACHER, Applications of the Theory of Matrices, Math. Comp. Review 43, 

v. 14, 1960, p. 284-285. 

This book is a translation and revision of the second volume of Gantmacher’s 
Theory of Matrices that was carried out by three people; namely, J. L. Brenner 
(named as the sole translator in the review under discussion), Mr. 8S. Evanusa and 
Prof. D. W. Bushaw. 


Muruan §8. Corrinetron, “Applications of the complex exponential integral,” 
Math. Comp., v. 15, 1961, p. 1-6. 


On p. 2, eq. (llc) should read Si(—x — iy) = ~—Si(x + iy) in place of 
Si(—x — ty) = Si(x + ty). 
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